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ABSTRACT 

Context. Density and velocity fluctuations on virtually all scales observed with modern telescopes show that molecular clouds (MCs) 
are turbulent. The forcing and structural characteristics of this turbulence are, however, still poorly understood. 
Aims. To shed light on this subject, we study two limiting cases of turbulence forcing in numerical experiments: solenoidal 
(divergence-free) forcing and compressive (curl-free) forcing, and compare our results to observations. 

Methods. We solve the equations of hydrodynamics on grids with up to 1024 3 cells for purely solenoidal and purely compressive 
forcing. Eleven lower-resolution models with different forcing mixtures are also analysed. 

Results. Using Fourier spectra and A-variance, we find velocity dispersion-size relations consistent with observations and indepen- 
dent numerical simulations, irrespective of the type of forcing. However, compressive forcing yields stronger compression at the same 
RMS Mach number than solenoidal forcing, resulting in a three times larger standard deviation of volumetric and column density 
probability distributions (PDFs). We compare our results to different characterisations of several observed regions, and find evidence 
of different forcing functions. Column density PDFs in the Perseus MC suggest the presence of a mainly compressive forcing agent 
within a shell, driven by a massive star. Although the PDFs are close to log-normal, they have non-Gaussian skewness and kurtosis 
caused by intermittency. Centroid velocity increments measured in the Polaris Flare on intermediate scales agree with solenoidal 
forcing on that scale. However, A-variance analysis of the column density in the Polaris Flare suggests that turbulence is driven on 
large scales, with a significant compressive component on the forcing scale. This indicates that, although likely driven with mostly 
compressive modes on large scales, turbulence can behave like solenoidal turbulence on smaller scales. Principal component analysis 
of G216-2.5 and most of the Rosette MC agree with solenoidal forcing, but the interior of an ionised shell within the Rosette MC 
displays clear signatures of compressive forcing. 

Conclusions. The strong dependence of the density PDF on the type of forcing must be taken into account in any theory using the 
PDF to predict properties of star formation. We supply a quantitative description of this dependence. We find that different observed 
regions show evidence of different mixtures of compressive and solenoidal forcing, with more compressive forcing occurring primar- 
ily in swept-up shells. Finally, we emphasise the role of the sonic scale for protostellar core formation, because core formation close 
to the sonic scale would naturally explain the observed subsonic velocity dispersions of protostellar cores. 

Key words, hydrodynamics — ISM: clouds — ISM: kinematics and dynamics — methods: numerical — methods: statistical — 
turbulence 



1. Introduction 

Studying the density and velocity distributions of interstellar 
gas provides essential information about virtually all physi- 
cal processes relevant to the dynamical evolution of the in- 
terstellar medium (ISM). Along with gravity, magnetic fields 
and the thermodynamics of the gas, supersonic turbulence 
plays a fundamental role in dete rmining the dens ity and ve- 
locity statistics of the ISM (e.g.. IScalo et al.lll998h . Thus, su- 
personic turbulence is considered a key process for star for- 



mation (Mac Low & Klessen 2004; Elmegreen & Scalo 2004; 
IScalo & Elmegreenll2004l;lMcKee & Ostrikerll2007l) . 

In this paper, we continue our analysis of the den- 
sity probability distribution function (PDF) obtained in nu- 
merical experiments of driven supersonic isothermal turbu- 
lence. Understanding the density PDF and its turbulent ori- 
gin is essential, because it is a key ingredient for an- 
alytical models of star formation: The turbulent density 

PDF is used to expla in the stellar initial mass func- 

tion dPadoan & Nordlundl l2002t iHennebelle & Chabrier| [2008. 
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20091). the star formation rate dKrumholz & McKed 120051: 
Krumholz et ail 120091: iPadoan & Nordlundll2009t). fliestar for- 



mation efficiency (Elmegreen 2008), and the Kennicutt- Schmidt 
relation on gala ctic scales dElmegreenll2002t iKravtsovl 12003: 
lTassisll2007b . In IFederrath et all d2008bl) . we found that super- 
sonic turbulence driven by a purely compressive (curl-free) force 
field yields a density PDF with roughly three times larger stan- 
dard deviation compared to solenoidal (divergence-free) turbu- 
lence forcing, which strongly affects the results obtained in these 
analytical models. Here, we want to compare our results for 
the density PDF to observations of column density PDFs (e.g., 
lGoodmanetal]|2009t). 

Moreover, in Federrath et al ] (I2009D we investigated the frac- 
tal density distribution of our two models with solenoidal and 
compressive turbulence forcing, which showed that compressive 
forcing yields a significantly lower fractal dimension (Df m 2.3) 
compared to solenoidal forcing (Df w 2.6). In the present contri- 
bution, we consider the scaling of centroid velocity increments 
computed for these models, and we compare them to observa- 
tions of the Polaris Flare by iHilv-Blant et al.l (120081) . We addi- 
tionally used principal component analysis and compared our 
results to observatio ns of the G2 16-2.5 (Maddalena's Cloud) and 
the Rosette MC bv lHever et ail d2006l) . 

Our results indicate that interstellar turbulence is driven by 
mixtures of solenoidal and compressive forcing. The ratio be- 
tween solenoidal and compressive modes of the turbulence forc- 
ing may vary strongly across different regions of the ISM. This 
provides an explanation for the apparent lack of correlation be- 
tween turbulent density and velocity dispersions found i n ob- 
servations (e.g.. iGoodman et al.l 12009; Pin eda et al .1120081) . We 
conclude that solenoidal forcing is more likely to be realised 
in quiescent regions with low star formation activity as in the 
Polaris Flare and in Maddalena's Cloud. On the other hand, in 
regions of enhanced stellar feedback, compressive forcing leads 
to larger standard deviations of the density PDFs, as seen in 
one of the subregions of the Perseus MC surrounding a central 
B star. Moreover, compressive forcing exhibits a higher scaling 
exponent of principal component analysis than solenoidal forc- 
ing. This higher scaling exponent is consistent with the mea- 
sured scaling exponent for the interior of an ionising shell in the 
Rosette MC. 

In §[2] we explain the numerical setup and turbulence forcing 
used for the present study. We discuss our results obtained using 
PDFs, centroid velocity increments, principal component analy- 
sis, Fourier spectrum functions, and A-variance analyses in § [3] 
mE]|6] and[7| respectively. In each of these sections, we compare 
the turbulence statistics obtained for solenoidal and compressive 
forcing with observational data available in the literature. In §|8] 
we discuss the possibility that transonic pre-stellar cores typi- 
cally form close to the sonic scale in a globally supersonic, tur- 
bulent medium. Section[9]provides a list of the limitations in our 
comparison of numerical simulations with observations. A sum- 
mary of our results and conclusions is given in § [10] 

2. Simulations and methods 



The piecewise parabolic method dColella & Woodward! 19841). 
implemented in the astro physical code FLASH3 dFrvxell et all 
2000; iDubev et al] [2008) was used to integrate the equations 
of hydrodynamics on three-dimensional (3D) periodic uniform 
grids with 256 3 , 512 3 , and 1024 3 grid points. Since isothermal 
gas is assumed throughout this study, it is convenient to define 



as the natural logarithm of the density divided by the mean den- 
sity (p) in the system. For isothermal gas, the pressure, P = pel, 
is proportional to the density p with the constant sound speed c s . 
The equations of hydrodynamics solved here are consequently 
given by 



ds 

— + (v • V)s = -V ■ v 
at 

^+(y. V)V = ~C?Vi+/, 



(2) 
(3) 



where v denotes the velocity of the gas. An energy equation 
is not needed, because the gas is isothermal. The assumption 
of isothermal gas is very crude, but may still provide an ad- 
equate physical appr oximation to the real thermodynamics i n 
dense molecular gas dWolfire et al.lll995l : lPavlovski et al.l l2006). 
We discuss further limitations of our simulations in § [9] The 
stochastic forcing term / is used to drive turbulent motions. 

2.1. Forcing module 

Equations (0 and (0 have been solved before in the con- 
text of molecular cloud dynamics, studying compressible tur- 
bulence with either solenoidal (divergence-free) forcing or 
with a 2:1 mixture of solen oidal to compressiv e modes in 
the turbulence forcing (e.g. , IPadoan et al 



1998; iMac Low et al 
2000: 



1997; Stone et al 



1998; Mac Low 1999; Klessen et al 



z.uuu[. iHe itsch et al. 2001; Klesse n 1200 It I Boldvrev et al. 
2002[ iLi et alj 120031; IPadoan et all 12004 Uappsen et all 120051; 
Ballesteros-Paredes et al.l 120061: iKritsuk et al.l 120071; iDib et al 
2008; Kissma nn et alJ 120081: lOffner et al.l 120081: ISchmidt eta! 



s = \n — 

(P) 



(1) 



2009). The case of a 2:1 mixture of solenoidal to compres- 
sive modes is the natural result obtained for 3D forcing, if 
no Helmholtz decomposition (see below) is performed. Then, 
the solenoidal modes occupy two of the three available spa- 
tial dimension s on average, while the c ompressive modes onl y 
occupy one (Elmegreen & Scaioll2004l: [Federrath et al. 2008b). 
In the present study, the solenoidal forcing case is thus also 
used as a control run for comparison with previous studies us- 
ing solenoidal forcing. However, we additionally applied purely 
compressive (curl-free) forcing and analysed the resulting turbu- 
lence statistics in detail. Each simulation at a resolution of 1024 3 
grid cells consumed roughly 100 000 CPU hours. Therefore, we 
concentrated on two extreme cases of turbulence forcing with 
high resolution: (1) the widely adopted purely solenoidal forc- 
ing (V ■ / = 0), and (2) purely compressive forcing (V x / = 0). 
However, we also studied eleven simulations at numerical res- 
olution of 256 3 in which we smoothly varied the forcing from 
purely solenoidal to purely compressive by producing eleven dif- 
ferent forcing mixtures. 

The forcing term / is often modelled with a spatially static 
pattern, for which the am plitude is adjus ted in time fol lowing the 
metho ds introduced by Mac Low et al.l dl998l) and IStone et all 
(fl998l) . This results in a roughly constant energy input on large 
scales. Other studies model the random f orcing term / such 
that it can vary in time and space (e.g . . IPadoan et al.1 [2004: 
IKritsuk et al.ll2007t IFederrath et dJl200lSMlScruriidt et alj|2 009). 
Here, we used the Ornstein-Uhlenbeck (OU) process to model 
/, which belongs to the latter type. The OU process is a well- 
defined stochastic process with a finite autocorrelation timescale. 
It can be used to excite turbul ent motions in 3D, 2D , and 
ID simulations as explained in lEswaran & Popd dl988l) and 
ISchmidt et all ([2006). Using an OU process enables us to control 
the autocorrelation timescale T of the forcing. The concept of 
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using the OU process to excite turbulence and the projections in 
Fourier space necessary to get solenoidal and compressive force 
fields are described below. 

The OU process is a stochastic differential equation describ- 
ing the evolution of the forcing term / in Fourier space (k- 
space): 

d/Gfc, f) = /«,(*) V}{k) AW{t) - f(k, t) j . (4) 

The first term on the right hand side is a diffusion term. This term 
is modelled using a Wiener process W(t), which adds a Gaussian 
random increment to the vector field given in the previous time 
step At. Wiener processes are random processes, such that 

W(i)-W(t-di) = N(0,df), (5) 

where A^O, df) denotes the 3D, 2D, or ID version of a Gaussian 
distribution with zero mean and standard deviation df. This is 
followed by a projection with the projection tensor 'Pj(k) in 
Fourier space. In index notation, the projection operator reads 

T\ } (Jfc) = ? Pjj (*) + (1 - P\j (*) = ( Sij + (1 - 20 , (6) 

where 6jj is the Kronecker symbol, and ff. = Sij - k/kj/k 2 and 
V\. - kikj/k 1 are the fully solenoidal and the fully compressive 
projection operators, respectively. The projection operator serves 
to construct a purely solenoidal force field by setting £ = 1 . For 
£ = 0, a purely compressive force field is obtained. Any combi- 
nation of solenoidal and compressive modes can be constructed 
by choosing f e [0, 1]. By changing the parameter we can 
thus set the power of compressive modes with respect to the total 
power of the forcing. The analytical ratio of compressive power 
to total power can be derived from equation (O by evaluating the 
norm of the compressive component of the projection tensor, 

|(W)?l/ = (W) 2 , (7) 

and by evaluating the norm of the full projection tensor 

|!Pf| 2 = 1 - 2( + . (8) 

The result of the last equation depends on the dimensionality 
D = 1,2,3 of the forcing, because the norm of the Kronecker 
symbol |<5y| = 1, 2 and 3 in one, two and three dimensions, re- 
spectively. The ratio of equations (0 and ([8]) gives the ratio of 
compressive forcing power Fi ong to the total forcing power F tot 
as a function of the parameter f : 

= (W) 2 

F tot i - 2<r + ■ 

Figure[T]provides a graphical representation of this ratio for 
the ID, 2D, and 3D case. For comparison, we plot numerical 
values of the forcing ratio obtained in eleven 3D and 2D hydro- 
dynamical runs with resolutions of 256 3 and 1024 2 grid points, 
in which we have varied the forcing parameter f from purely 
compressive forcing (£ = 0) to purely solenoidal forcing (f = 1) 
in the range £ = [0,1], separated by A£ = 0.1. Note that a natural 
mixture of forcing modes is obtained for ( = 0.5, which leads 
to Fi ong /F tot = 1/3 for 3D turbulence, and Fi ong /F tot = 1/2 for 
2D turbulence. A simple way to understand this natural ratio is 
to consider longitudinal and transverse waves. In 3D, the longi- 
tudinal waves occupy one of the three spatial dimensions, while 




0.0 0.2 0.4 0.6 0.8 1.0 



Fig. 1. Ratio of compressive power to the total power in the 
turbulence force field. The solid lines labelled with ID, 2D, 
and 3D show the analytical expectation for this ratio, equa- 
tion (0, as a function of the forcing parameter £ for one-, two- 
and three-dimensional forcing, respectively. The diamonds and 
squares show results of numerical simulations in 3D and 2D 
with £ = [0,1], separated by A£ = 0.1. Those models were 
run at a numerical resolution of 256 3 and 1024 2 grid points 
in 3D and 2D, respectively. The two extreme forcing cases of 
purely solenoidal forcing = 1) and purely compressive forc- 
ing = 0) are indicated as "sol" and "comp", respectively. Note 
that in any ID model, all power is in the compressive component, 
and thus F\ ong /F toi = 1, independent of 



the transverse waves occupy two of the three on average. Thus, 
the longitudinal (compressive) part has a power of 1/3, while the 
transverse (solenoidal) part has a power of 2/3 in 3D. In 2D, the 
natural ratio is 1/2, because longitudinal and transverse waves 
are evenly distributed in two dimensions. 

The second term on the right hand side of equation (O is 
a drift term, which models the exponentially decaying corre- 
lation of the force field with itself. Thus, the autocorrelation 
timescale of the forcing is denoted by T. We set the autocor- 
relation timescale equal to the dynamical timescale T = L/(2V) 
on the scale of energy injection, where L is the size of the com- 
putational domain, V - c s M and M ~ 5.5 is the RMS Mach 
number in all runs. The autocorrelation timescale is therefore 
equal to the decay time constant in supersonic hydrodynamic 
and magnetohydrodynamic turbu lence driven on large scales 
dStone et alJl998HMac Lowlll999h . The forcing amplitude f (k) 
is a paraboloid in 3D Fourier space, only containing power on 
the largest scales in a small interval of wavenumbers 1 < \k\ < 3 
peaking at k = 2, which corresponds to half of the box size L/2. 
Th e effects of varying the scale of energy i n put were investigated 
bv lMac Low! d!999l).lKlessen et alj d2000l) . Heitsch etafl d2001l) 
and IVazquez-Semadeni et alJ (l2003h . Here, we consider large- 
scale stochastic forcing, which is closer to the observa tional data 
(e.g., Ossenkopf & Mac Lowll2002l:lBrunt et alj2009l) . This type 
of forcing models the kinetic energy input from large-scale tur- 
bulent fluctuations, breaking up into smaller structures. Kinetic 
energy cascades down to smaller and smaller scales, and thereby 
effectively drives turbulent fluctuations on scales smaller than 
the turbulence injection scale. 

We have verified that our results are not sensitive to the gen- 
eral approach of using an Ornstein-Uhlenbeck process for the 
turbulence forcing. For instance, we have used an almost static 
forcing pattern, which is obtained in the limit T — > oo in test sim- 
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ulations. We have furthermore checked that the particular choice 
of Fourier amplitudes did not affect our results by using a band 
spectrum instead of a parabolic forcing spectrum. Varying these 
parameters did not strongly affect our results. In contrast, chang- 
ing f from f = 1 (solenoidal forcing) to £ = (compressive forc- 
ing) always led to significant changes in the turbulence statistics. 

2.2. Initial conditions and post-processing 

Starting from a uniform density distribution and zero veloci- 
ties, the forcing excites turbulent motions. Equations © and Q 
have been evolved for ten dynamical times T, which allows us 
to study a large sample of realisations of the turbulent flow. 
Compressi ble turbulence reached a statistically invariant state 
within 2 T ( Federrath et al. 2009). This allows us to average all 
statistical measures over 8 T separated by 0. 1 T in the fully de- 
veloped regime. We are thus able to average over 81 different 
realisations of the turbulence to improve statistical significance. 
The 1-cr temporal fluctuations obtained from this averaging pro- 
cedure are indicated as error bars for the PDFs, centroid velocity 
increments, principal component analysis, Fourier spectra and 
A-variance analyses in the following sections and in all figures 
showing error bars throughout this study. The forcing amplitude 
was adjusted to excite a turbulent flow with an RMS Mach num- 
ber M ~ 5.5 in all cases. We use the RMS Mach number as 
the control parameter, because this dimensionless number deter- 
mines most of the physical properties of scale-invariant turbulent 
flows and is often used to derive important flow statistics such as 
the standard deviation of the density distribution. However, in 
the next section we show that the latter depends sensitively on 
the turbulence forcing parameter f as well. 

Figure [2] (top panels) shows column density fields projected 
along the z-axis from a randomly selected snapshot at time 
t = 2 T in the regime of fully developed, statistically station- 
ary turbulence for solenoidal (left) versus compressive forcing 
(right). This regime was reached after 2 dynamical times T, 
which is shown in Figure[3]for the minimum and maximum log- 
arithmic densities s (top panel) and RMS curl and divergence of 
the velocity field (bottom panel) as a function of the dynamical 
time. It is evident that compressive forcing produces higher den- 
sity contrasts, resulting in higher density peaks and bigger voids 
compared to solenoidal forcing. 



3. The probability density function of the gas 
density 

It is interesting to study the probability distribution of tur- 
bulent density fluctuations, because it is a key ingredient 
for many analytical models of star form ation: it is used to 
expla i n the stellar initial mass func t ion (jPadoan & Nordlundl 
2002] iHennebelle & Cha brier 120081 l2009ir the star forma- 



tion rate dKrumholz & McKed l2005t IKrumholz etail l2009t 
Padoan & Nordlundl 12009), the star formation efficiency 
(Elmegreen 2008), and the Kennicutt-Schmidt relation on galac- 
tic scales (Elmegre enl2002HKravtsovll2003HTassis1l2007l) . 

The probability to find a volume with gas density in the 
range [p, p + dp] is given by the integral over the volume- 
weighted probability density function (PDF) of the gas density: 
jP +d P p p (p'~)dp' Thus, the PDF p describes a probability den- 
sity, which has dimensions of probability divided by gas density 
in the case of p p (p). By the same definition, p s (s) denotes the 
PDF of the logarithmic density s = ln(p/ (p)). 



Figure [4] presents the comparison of the time-averaged 
volume-weighted density PDFs p s (s) obtained for solenoidal 
and compressive forcing. The linear plot of p s (s) (top panel) dis- 
plays the peak best, whereas the logarithmic representation (bot- 
tom panel) reveals the low- and high-density wings of the dis- 
tributions. Three different fits to analytic expressions (discussed 
below) are shown as well. 

3. 1 . The density PDF for solenoidal forcing 

In numerical experiments of driven supersonic isother- 
mal turbulence with solenoidal and/or wea k ly compres- 

[1994; Padoa n et alJ 
Nordlund & Padoan 
IPadoan etalJl2004t 



sive forcing (e.g.. IVaz quez-Semadeni 



1 1997b IStoneetaU Il998t iMac Low! 199S 
I 1999L iBoldvrev et al] l2002t iLi et alJ l2003t 
iKritsuk et al]|2007t iBeetz et a l. 2008), but also in decaying tur- 
bulen ce (e.g.. lOstriker et alJll999t lKlessenll2000l: lOstriker et all 
1200 U lCH over & Mac Lo w 2007b) it was shown that the density 
PDF p s is close to a log-normal distribution, 



Ps ds 



exp 



(s-(s)) 2 



2cr 2 



ds 



(10) 



where the mean (s) is related to the standard deviation cr s by 
(s) = —cr 2 J2 due to the constraint of mass conservation (e.g., 
IVazquez- Semaden il 1 9941) : 



exp (s) p s ds = I p p p dp = (p) 

oo Jo 



(11) 



Equation dTTJ simply states that the mean density has to be re- 
covered. This constraint together with the PDF normalisation, 



p s ds= I 

oo Jo 



Pp dp = 1 



(12) 



must always be fulfilled for any density PDF whether log-normal 
or non-Gaussian. 

From our simulations, we obtain density PDFs in agreement 
with log-normal distributions for solenoidal forcing. The log- 
normal fit using equation ( TTOb is shown in Figure |4] as dashed 
lines. However, the PDF is not perfectly log-no rmal, i.e., there 
are w eak non-Gaussian contributions (see also, iDubinski et al .1 
1 19951) . especially affecting the wings of the distribution. The 
strength of these non-Gaussian features is quantified by com- 
puting higher-order moments (skewness and kurtosis) of the dis- 
tributions. The fir st four standardised central moments (see, e.g., 
IPress et aT1ll986l) of a discrete dataset {q} with elements are 
defined as 



1 N 

mean : (q) = — ^ q t 

i=\ 

dispersion: a q = yj{(q - (q)) 2 ) 

((q-{q)) 3 ) 



(13) 



skewness : S a 



((q-(q)) 4 ) 



kurtosis : 7C = 

q a 4 

Note that in our definition of the kurtosis (also called flatness), 
the Gaussian distribution has 7C = 3. We have computed the first 
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Fig. 2. Maps showing density (top), vorticity (middle) and divergence (bottom) in projection along the z-axis at time t — 2 T as 
an example for the regime of statistically fully developed, compressible turbulence for solenoidal forcing (left) and compressive 
forcing (right). Top panels: Column density fields in units of the mean column density. Both maps show three orders of magnitude 
in column density with the same scaling and magnitudes for direct comparison. Middle panels: Projections of the modulus of 
the vorticity |V x v|. Regions of intense vorticity appear to be elongated filamentary structures often coinciding with positions of 
intersecting shocks. Bottom panels: Projections of the divergence of the velocity field V ■ v showing the positions of shocks. Negative 
divergence corresponds to compression, while positive divergence corresponds to rarefaction. 
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t / T 



1 024 3 sol 
1 024 3 comp 
log-normal fit 
skew — log — normal fit 
skew-kurt-log-normal fit 



volumetric density PDFs 




Fig. 3. Top panel: Minimum and maximum logarithmic den- 
sity s — In (pi {p}) as a function of the dynamical time T. 
Note that compressive forcing yields much stronger compres- 
sion and rarefaction compared to solenoidal forcing, although 
the RMS Mach number is roughly the same in both cases 
(see|Federrath et al.ll2009l Fig. 2). Bottom panel: RMS vorticity 
<(Vxv) 2 ) 1/2 and RMS divergence <(V-v) 2 } 1/2 as a function of the 
dynamical time. Within the first 2 T, a statistically steady state 
was reached for both solenoidal (sol) and compressive (comp) 
forcing. This allows us to average statistical measures (prob- 
ability density functions, centroid velocity increments, princi- 
pal component analysis, Fourier spectra and A-variances) in the 
range 2 < t/T < 10 to improve statistical significance of our 
results and to estimate the amplitude of temporal fluctuations 
(snapshot-to-snapshot variations) between different realisations 
of the turbulence. 



four statistical moments of the volumetric density PDFs shown 
in Figure |U The results are summarised in Table Q] The 1-cr er- 
ror given for each statistical moment was obtained by averaging 
over 81 realisations of the turbulence as described in § 12.21 Both 
solenoidal and compressive forcing yield density PDFs with de- 
viations from the Gaussian 3rd order (skewness S = 0) and 4th 
order (kurtosis "7C = 3) moments. 

3.2. The density PDF for compressive forcing 

Contrary to the solenoidal case, the PDF obtained for compres- 
sive forcing is not at all well fitted with the perfect log-normal 
functional form (dashed line in Figure |4] for compressive forc- 
ing). Due to the constraints of mass conservation (eq. [TTI) and 
normalisation (eq. fT2l) . the peak position and its amplitude can- 
not be reproduced simultaneously. The skewness and kurtosis 
for the compressive forcing case are also listed in Table [1] Non- 
Gaussian values of skewness and kurtosis, i.e., higher-order mo- 
ments require modifications to the analytic expression of the log- 
normal PDF given by equation ([Tol l. A first step of modification 



Fig. 4. Volume-weighted density PDFs p(s) of the logarithmic 
density s = ln(p/ (p)) in linear scaling (top panel), which dis- 
plays the peak best, and in logarithmic scaling (bottom panel) to 
depict the low- and high-density wings. The PDF obtained from 
compressive forcing (1024 3 comp) is significantly wider than the 
solenoidal one (1024 3 sol). The peak is shifted to lower values 
of the logarithmic density s, because of mass conservation, de- 
fined in eq. (fTTT i. The density PDF from solenoidal forcing is 
compatible with a Gaussian distribution. However, there are also 
non-Gaussian features present, which are associated with inter- 
mittency effects. These are more prominent in the density PDF 
obtained from compressive forcing, exhibiting statistically sig- 
nificant deviations from a perfect log-normal (fit using eq. [10] 
shown as dashed lines). A skewed log-normal fit (dash-dotted 
lines) given by eq. (fT4l provides a better representation, but still 
does not fit the high-density tail of the PDF obtained for com- 
pressive forcing. Both the PDF data obtained from solenoidal 
and compressive forcing are best described as log-normal distri- 
butions with higher-order corrections defined in eq. ( fTTl l. which 
take into account both the non-Gaussian skewness and kurtosis 
of the distributions. These fits are shown as solid lines (skew- 
kurt-log-normal fit). The first four standardised moments defined 
in equations ( fT3l of the distributions in p and s are summarised 
in Table Q] together with the fit parameters. The grey shaded re- 
gions indicate 1-cr error bars due to temporal fluctuations of the 
distributions in the regime of fully developed, supersonic turbu- 
lence. A total number of 1024 3 x81 ~ 10" data points contribute 
to each PDF. 



is to allow for a finite sk ewness, which is possible with a skewed 
log-normal distribution riAzzaiinill 19851) 



P(s) 



71 U) 



exp 



(s-tY 



2(D 2 



X(s-f)a-/w / f 2 \ 

. exp b) d? 



(14) 
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Table 1. Statistical moments and fit parameters of the PDFs of the volumetric density p for solenoidal and compressive forcing 
shown in Fig. |4] 



Solenoidal Forcing Compressive Forcing 



Standardised Moments 


(o) 


1.00 ± 


0.00 


1.00 


± 0.00 


of p 




1.89 ± 


0.09 


5.86 


± 0.96 






9.03 ± 


1.06 


26.7 


± 10.1 




'K 


211. ± 


69.8 


1720 


± 2000 


Standardised Moments 


(s) 


-0.83 ± 


0.05 


-3.40 


± 0.43 


of v - In (nl (oW 


u S 


1.32 ± 


0.06 


3.04 


± 0.24 




s s 


-0.10 ± 


0.11 


-0.26 


+ 0.20 






3.03 ± 


0.17 


2.91 


+ 0.43 


Skewed Log-normal Approximation 




0.010 ± 


0.050 


-0.048 


±0.133 


using equation d!4t 


CO 


1.562 ± 


0.035 


4.712 


±0.193 




a 


-0.911 ± 


0.064 


-2.163 


±0.173 


4th Order Approximation 


"0 


- 1.3664 ± 


0.0091 


-2.5014 


± 0.0259 


(including skewness and kurtosis) 


fill 


-0.4592 ± 


0.0064 


-0.3437 


±0.0132 


using equation dl7t 


a 2 


-0.3067 ± 


0.0052 


-0.0831 


± 0.0030 




a 3 


-0.0073 ± 


0.0011 


-0.0065 


±0.0011 




a \ 


-0.0002 ± 


0.0005 


-0.0004 


± 0.0001 



where a, £ and u> are fit parameters. Defining 6 — a/ Vl + a 2 , 
the first four standardised central moments of the distribution are 
linked to the parameters a, £ and to, such that 



mean : 



(s) = £ + co6^2/tt 



vl/2 



dispersion: cr s = co(l-26 2 /jrj 

4-71 (6y/2/n) 3 



(15) 



skewness 



kurtosis 



2 (1 -2£ 2 /V> 3/2 

2(n - 3)(6 V27i) 4 
(1 - 2d 2 In) 2 ' 



Skewed log-normal fits are added to Figure [4] as dash-dotted 
lines and the corresponding fit parameters are given in Table Q] 
However, for a skewed log-normal distribution, the kurtosis is 
a function of the skewness, since the skewness and kurtosis in 
equations (fT3T > both depend on the same parameter 6 only. 

Better agreement between an analytic functional form and 
the measured PDF can be obtained, if the actual kurtosis of the 
data is taken into account as an independent parameter in the 
analytical approach. The fundamental derivation of a standard 
Gaussian distribution is given by 



In p(s) — oo + CL\S + ci2S 



(16) 



where one parameter is constrained by the normalisation and the 
two remaining ones are determined by the mean and the disper- 
sion. We can extend this to a modified Gaussian-like distribution 
by including higher-order moments: 

p(s) = exp [ao + ci\s + 02s 2 + 03s 3 + 04s 4 + (5(s 5 )] . (17) 

Here, the expansion is stopped at the 4th moment. One parameter 
is again given by the normalisation, and the remaining four pa- 
rameters are related to the mean, dispersion, skewness and kur- 
tosis. Fits obtained with this formula are included in Figure [4] 
as solid lines. The fit parameters are listed in Table Q] This 



new functional form is in good agreement with the data from 
solenoidal and compressive forcing, fitting both the peak and the 
wings very well. They follow the constraints of mass conser- 
vation and normalisation given by equations (fTTb and (1121) . We 
have computed the first four moments of the fitted function and 
find very good agreement with the first four moments of the ac- 
tual PDFs. 

The fitted parameters 03 and 04, which represent the higher- 
order terms tend to zero compared to the standard Gaussian pa- 
rameters at), a\ and «2 (see Table[TJ. This means that the higher- 
order corrections to the standard Gaussian are small. However, 
we point out that they are absolutely necessary to obtain a good 
analytic representation of the PDF data, given the fact that equa- 
tions (fTTT i and ( TTZl i must always be fulfilled and that the analytic 
PDF should return the correct values of the numerically com- 
puted moments of the measured distributions. 

In the various independent numerical simulations mentioned 
above, the density PDFs were close to log-normal distributions 
as in our solenoidal and compressive forcing cases. However, 
most of these studies also report considerable deviations from 
Gaussian PDFs, which affected mainly the low- and high-density 
wings of their distributions. These deviations can be associated 
with rare events caused by strong intermittent fluctuations during 
head-on collisions of strong shock s and oscillations in very low- 
densi t y rarefaction wave s (e.g., iPassot & Vazquez-Semadenil 
1998: iRritsuk et al1l2007h . The pronounced deviations from the 
log-normal shape of the density P DF for compre s sively driven 
turbulence were also discussed by Schmidt et al] (120091) . Even 
stronger deviations from log-normal PDFs we re reported in 
strongly self-gravitating turbulent systems (e.g ., iKlessenfe OOO; 
iFederrath et al.ll2008at Kainulainen eTaT]l2009h . 

Intermittency is furthermore inferred from obser- 
vations, affecting the wi ngs of molecular line profiles 
dFalgarone & Phillipsl 1 19901). and the stati sti cs of centroid 
veloci ty increments dHilv-Blant et al.l 120081) . iGoodman et al.l 
(2009) measured column density PDFs using dust extinction 
and emission, as well as molecula r lines of gas in t he Pe rseus 
MC. Using dust extinction maps, Lombardi et al. (2006) ob- 
tained the column density PDF for the Pipe nebula. The PDFs 
found in these studies roughly follow log-normal distributions. 
However, deviations from perfect log-normal distributions are 
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clearly present in the density PDFs obtained in these studies. 
They typically exhibit non-Gaussian features. For instance, 
lLombardi et al.l (12006) had to apply combinations of multiple 
Gaussian distributions to obtain good agreement with the 
measured PDF data. 



3.3. Density-Mach number correlation and signatures of 
intermittency in the density PDFs 

As discussed by iPassot & Vazquez-Semadenil d 1998b . a 
Gaussian distribution in the logarithm of the density, i.e., 
a log-normal distribution in p is expected for supersonic, 
isothermal turbulent flows. The fundamental assumption be- 
hind this model is that density fluctuations are built up by 
a hierarchical process. The local density p(r, t) at a given 
position r is determined by a Markov process, i.e., by the 
product p(t„) = S(t n -i)p(t n -i) = . . . = 6(t )p(t ) of a large 
number of independent ra ndom fluctuations 6{t n ) > in time 
dVazquez- Semadenil |1994|) . If these fluctuations were indeed 
independent, the quantity s = ln(p/ (p)) would be determined 
by the sum of this large number of local fluctuations and the dis- 
tribution in s becomes a Gaussian distribution according to the 
central limit theorem. Since the equations (0 and (0 are invari- 
ant under the transformation s — > s + sq for an arbitrary constant 
So, the random variable s(t n ) should be independent of the local 
Mach number, and independent of t he density at previous times 
?„_!, t n -2, ■■■,?()■ As pointed out b y Vazquez-Semadenil d 19941) 
and IPassot & Vazquez-Semadenil (1998), this independence 
breaks down in strong shocks and density extrema, because 
«o cannot be arbitrarily high due to mass conservation, and 
an upper boundary s + exists. In consequence, if s+ is reached 
locally, the density cannot increase anymore by a subsequent 
fluctuation, and the next density is not independent of the 
previous timestep, causing the fundamental assumption to break 
down. This also applies to strong rarefaction waves, because 
creating shocks always produces strongly rarefied regions 
outside the shock. 

When the fundamental assumption breaks down, den- 
sity and velocity statistics are expected to become correlated 
dVazquez-Semadenil [19941: IPassot & Vazquez-Semadenil ! 1998; 
Kritsu k et al.l l 20071) . Since in isothermal gas, the sound speed is 
constant, this translates directly into Mach number-density cor- 
relations. The average local Mach number M = v/c s may there- 
fore exhibit some dependence on the average local density. For 
instance, it is intuitively clear that head-on collisions of strong 
shocks produce very high density peaks. In the stagnation point 
of the flow, the local velocity and consequently the local Mach 
number will almost drop to zero. The time evolution of the max- 
imum and minimum de nsity in Figure[3]sho ws these intermitten t 
fluctuations (see also, lPorteretal.lll992bt iKritsuk et all 120071) . 
The intermittent phenomenon corresponds to the situation ex- 
plained above, for which s + might have been reached, and some 
dependence of the Mach number on density is expected. 

In real molecular clouds, the maximum densities are sim- 
ilarly bounded, and cannot reach infinitely high values, either. 
This is-unlike the finite resolution constraints in simulations- 
because the gas becomes optically thick at a certain density 



Table 2. Standard deviations of the density PDFs as a function 
of numerical resolution for solenoidal and compressive forcing 
shown in Fig. [6] 



(p > 10 14 g cm 3 ), an d cannot cool effici ently anymore (e.g., 
lLarsonlll969t lPenstod[l969T: iLarsonl 120051: iJappsen et alJl2005l 
and references therein). The gas is not close to isothermal any- 
more in this regime, and adiabatic compression induced by tur- 
bulent motions remain finite in real molecular clouds. Thus, the 
reason for the breakdown of the density-Mach number indepen- 
dence is different in simulations and observations, but it might 



Grid Res. 


Solenoidal Forcing 


Compressive Forcing 




cr p cr, 


Cp CTj 


256 3 ... 


1.79 ±0.08 1.36 ±0.07 


5.66 ±0.79 3.09 ±0.21 


512 3 ... 


1.89 ±0.10 1.35 ±0.05 


5.59 ±0.67 3.15 ±0.34 


1024 3 ... 


1.89 ±0.09 1.32 ±0.06 


5.86 ±0.96 3.04 ±0.24 



still be fundamental for the deviations from a log-normal PDF. 
Moreover, the existence of a characteristic scale may lead to a 
breakdown of the hierarchical model, and thus to a breakdown 
of the fundamental assumption. The scale at which supersonic 
turbulence becomes subsonic is such a scale. This scale is called 
the sonic scale, and is discussed later in §[8] 

We have computed the probability distributions for Mach 
number-density correlations. Figure [5] shows the volume- 
weighted correlation PDFs of local Mach number M ver- 
sus density p. Although the correlation between density and 
Mach number is weak as expected for isothermal turbulenc e 
(IVazquez-Semadeni|[T994l IPassot & Vazquez-Semadenii ri998). 
these two quantities are not entirely uncorrelated, which may ex- 
plain the deviations from perfect log-normal distributions. There 
is a weak trend for high-density regions to exhibit lower Mach 
numbers on average. Power-law estimates for densities above the 
mean logarithmic density indicate Mach number-density corre- 
lations of the form M(p) cc p-° M for solenoidal and M(p) oc 
p-0.05 f or com p ress i v e forcing. A similar power law exponent 
can be obtained from lKritsuk et al.1 (120071 Fig. 4). 



3.4. Numerical resolution dependence of the density PDFs 

The high-density tails of the PDFs in Figure |4] are not perfectly 
fit, even when the skewness and kurtosis are taken into account. 
This is partly due to non-zero 5th, 6th and higher-order mo- 
ments in the distributions, and partly because our numerical res- 
olution is insufficient to sample the high-density tail perfectly. 
Figure|6]shows that even at a numerical resolution of 1024 3 grid 
points, the high-density tails are not converged in both solenoidal 
and compressive forcing and tend to underestimate high den- 
sities. This limitation is shared among all turbulence simula- 
tions (see, for instan ce, the turbulence comparison project by 
Kitsionas et "all l2009l) . since the strongest and most intermittent 
fluctuations building up in the tails will always b e truncated due 
to limited numerical resolution (see also [Hennebelle & Audit! 
120071; iKowal et al.ll2007tlPrice & Federrathl2010l) . However, the 
peak and the standard deviation of the PDFs are reproduced quite 
accurately at a resolution of 256 3 . Table [2] shows the values of 
the linear standard deviation <x p and logarithmic standard de- 
viation <t s for numerical resolutions of 256 3 , 512 3 and 1024 3 . 
There appears to be no strong systematic dependence of the 
standard deviations on the numerical resolution for resolutions 
above 256 3 . The statistical fluctuations are the dominant source 
of uncertainty in the derived values of the standard deviations. 
It should be noted however that we have tested only the case 
of an RMS Mach number of about 5-6 here. There might be a 
stronger resolution dependence for higher Mach numbers, due 
to the stronger shocks produced in higher Mach number turbu- 
lence, which should be tested in a separate study. 
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Fig. 5. Volume-weighted correlation PDFs of local Mach number M versus logarithmic density s for solenoidal (left) and compres- 
sive forcing (right). Adjacent contour levels are spaced by 0.25 dex in probability density. Density and Mach number exhibit a very 
weak, but non- zero correlation in both forcing cases, which provides an explanation for the non -Gaussian features in the density 
PDFs of Fig.l4l(Vazquez-Semade mTl994t IPassot & Vazquez-Semadeni 1998: lKritsuk et aLll2007l) . The two solid lines, intersecting 
the maxima of both distributions show the mean Mach number as a function of the logarithmic density s = ln(p/ (p)). The tendency 
for high-density gas having lower Mach numbers on average is indicated as power laws in the high-density parts of the distributions. 
This suggests that the Mach number M(p) oc p -0 06 for solenoidal and M(p) cc p-° 05 for compressive forcing. 



comp 




1024 3 


comp 


512 3 


comp 


256 3 


comp 


1024 3 


sol 


512 3 


sol / 


256 3 


sol 



Fig. 6. Density PDFs at numerical resolutions of 256 3 , 512 3 and 
1024 3 grid cells. The PDFs show very good overall convergence, 
especially around the peaks. Table [2] shows that the standard 
deviations are converged with numerical resolution. The high- 
density tails, however, are not converged even at a numerical 
resolution of 1024 3 grid points, indicating a systematic shift to 
higher densities with resolution. Thi s limitation is shared among 
all turbulence simulations (see also, [H ennebelle & AuditH2007t 
iKitsionas et ai]|2009t IPrice & Federrafhll2010l) . The low-density 
wings are subject to strong temporal fluctuations due to inter- 
mittent bursts caused by head-o n collisions of shock s followed 
by strong rarefaction waves (e.g jKritsuk et alj2007l) . The inter- 
mittency causes deviations from a perfect Gaussian distribution 
and accounts for non-Gaussian higher-order moments (skewness 
and kurtosis) in the distributions. 



3.5. The column density PDFs and comparison with 
observations 

The strong difference between the statistics of the solenoidal and 
compressive forcing cases seen in the PDFs of the volumetric 
density shown in Figure |4]is reflected by the corresponding col- 



umn density PDFs. The time-averaged and projection-averaged 
column density PDFs are shown in Figure [7] Analogous to 
Table Q]for the volumetric density PDFs, we summarise the sta- 
tistical quantities and fit parameters for the column density PDFs 
in Table [3] The main results and conclusions obtained for the 
volumetric density distributions also hold for the column density 
distributions. Compressive forcing yields a column density stan- 
dard deviation roughly three times larger than solenoidal forc- 
ing. The relative difference between solenoidal and compres- 
sive forcing is thus roughly the same for the volumetric and 
the column density distributions. However, the absolute values 
are lower for the column density distributions compared to the 
volumetric density distributions. The reason for this is that by 
computing projections of the volumetric density fields, density 
fluctuations are effectively averaged out by integration along the 
line-of-sight, and as a consequence, the column density disper- 
sions become smaller compared to the corresponding volumetric 
density dispersions. 

The small inset in the upper right corner of Figure|7]addition- 
ally shows the column density PDFs computed along the z-axis 
at one single time t — 2 T corresponding to the map shown in 
Figure|2] This figure shows the effect of studying one realisation 
only, without time- and/or projection-averaging. This is interest- 
ing to consider, because observations can only measure column 
density distributions at one single time. Improving the statis- 
tical significance would only be possible by studying multiple 
fields and averaging in space r ather than in time invok ing the er- 
godic theorem as suggested by Goodman et al. (2009). However, 
even by studying one turbulent realisation only, the difference 
between solenoidal and compressive forcing is recovered from 
the dispersions of the distributions. We therefore expect that us- 
ing observations of column density PDFs, one can distinguish 
purely solenoidal from purely compressive forcing by measur- 
ing t he dispersion of the co lumn density PDF. 

iGoodman et al.l 12009) provided measurements of the col- 
umn density PDFs in the Perseus MC obtained with three dif- 
ferent methods: dust extinction, dust emission, and 13 CO gas 
emission. Although systematic differences were found between 
the three methods, they conclude that in general, the measured 
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Table 3. Same as Table[T] but for the PDFs of the column density X shown in Fig. [7] 



Solenoidal Forcing Compressive Forcing 



Standardised Moments 


(T> 


1.00 


± 0.00 


1.00 


± 0.00 


of £ 


(Tx 


0.47 


± 0.05 


1.74 


± 0.43 




5v 


1.38 


±0.38 


4.57 


± 1.44 




'A"v 


6.32 


±2.20 


36.8 


± 24.3 


Standardised ^Moments 


(n) 
VI/ 


-0.10 


± 0.02 


-1.00 


± 0.33 


of n = In (X/ (E» 

\ji if 111 V^ - / J 


" 1} 


0.46 


± 0.06 


1.51 


± 0.28 




s 


-0.04 


± 0.30 


-0.17 


± 0.29 




>I 


2.97 


±0.40 


2.69 


±0.45 


Skewed Log-normal Approximation 


€ 


0.180 


± 0.088 


0.717 


±0.102 


using equation d 1 4b 


CO 


0.539 


± 0.058 


2.392 


±0.160 




a 


-0.878 


± 0.342 


-2.371 


± 0.300 


4th Order Approximation 


«o 


-0.1524 


±0.0451 


-1.4547 


± 0.0532 


(including skewness and kurtosis) 


a\ 


-0.3900 


±0.1080 


-0.2902 


± 0.0417 


using equation d 1 71> 


"2 


-2.4643 


±0.1994 


-0.2669 


± 0.0259 






-0.1748 


±0.1469 


-0.0370 


± 0.0080 




a 4 


0.0204 


±0.1239 


-0.0035 


± 0.0024 



column density PDFs are close to, but not perfect log-normal 
distributions, which is consistent with our results. They further- 
more provided the column density PDFs and the column den- 
sity dispersions for six subregions in the Perseus MC. The dif- 
ference between the dispersions measured for these subregions 
is not as large as the difference between purely solenoidal and 
purely compressive forcing. The largest difference in the col- 
umn density dispersio ns among the six subregions found by 
iGoodman et al.1 (120091) is only about 50% relative to the aver- 
age column density dispersion measured in the Perseus MC. 
This indicates that both purely solenoidal and purely compres- 
sive forcing are very unlikely to occur in nature. On the other 
hand, a varying mixture of solenoidal and compressive modes 
close to the natural mixture of 2:1 can easily explain the 50% 
difference in density dispersion meas ured among the d ifferent 
regions. In particular, the Shell region (Ridge et alj2006l) . which 
surrounds the B star HD 278942 exhibits the largest density 
disper sion among all the subregions studied by IGoodman et al.l 
(2009), although its velocity dispersion is rather small compared 
to the others. This indicates that turbulent motions may be driven 
compressively rather t han solenoidally within the Shell region. 
IGoodman et al.l d2009l) indeed mentioned that the gas in the Shell 
is dominated by an "obvious driver", skewing the column den- 
sity distribution towards lower values compared to the other re- 
gions. Due to the constraints of mass conservation (eq. ITTb and 
normalisation (eq.[T2l. both the peak position and the peak value 
of the PDF skew to lower values, if the density dispersion in- 
creases (see Figure |7). Taken together, this suggests that the 
Shell in the Perseus MC represents an example of strongly com- 
pressive turbulence forcing rather than purely solenoidal forcing. 



3.6. The forcing dependence of the density dispersion-Mach 
number relation 

In iFederrath et al.l (12008b). we inve stigated the density 
dispersion-Mach number relat ion ([Padoan etaLl \\WR 
iPassot & Vazquez-Sem adeni 1998J) 1 , 



(Tp_ 



bM 



(18) 



This relation was also investigated in iKowal etail ([2007, 
Fig. 11), indicating that the standard deviation of turbulent 
density fluctuations, cr p is directly proportional to the sonic 
Mach number in the supersonic regime. It must be noted, 
however, that it was only d irectly tested for r ather low RMS 
Mach numbers, M < 2.5 dKowal et al.1 120071) and M < 3.5 
(Pass ot & Vazquez- Semadenil l 19981) . compared to typical Mach 
numbers in molecular clouds. If additionally a log-normal PDF, 
equation (TlOb is assumed, then equation ( fT8b can be expressed as 



=ln(l +b 2 M 2 ) 



(19) 



with th e same parameter b (Pad oan et alJll997T : Federrat h et al.l 
l2008bh . 

We begin our discussion of the forcing dependence of 
the density di spersion-Mach numbe r rel ation with a prob- 
lem rais ed by iMac Low et al.l d2005l) and iGlover & Mac Lowl 
(l2007bl) . lMac Low et alJ(l2005l) and lGlover & Mac Lowl(l2007bl) 
cla imed that the density dispersion-Mac h number relation found 
by IPassot & Vazquez-Semadenil (1 19981) . cr s = bM (which is a 
Taylor expansion of eq. [T9]for small RMS Mach numbers), with 
b as 1 did not at all fit their results for pressure and de nsity PDFs, 
while equation ( fT9T > with b « 0.5 (Pado an et al.l 19 97) provided a 
much better representation of their data. The main difference in 
the de nsity dispersion-Mach number relations b y IPadoan et al. 
d 19971) and IPassot & Vazquez-Semadenil dj998) i s the propor- 
tionali ty co nstant b. It is & 0.5 and b & 1 in iPadoan et al.l 
(1 19971) and IPassot & Vazquez-Semadenil d!998l) . respectively. 



Note that lPassot & Vazquez-Semadeni d!998l) suffers from a num- 
ber of typo graphical errors as a re sult of last-minute change of notation. 



Please see 
tions 



Mac Low et alj j2005i, footnote 5) for a number of correc- 
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4— 1024 3 sol 
x— 1024 3 comp 
— log-normal fit 
-- skew — log — normal fit 
— skew-kurt-log-normal fit 



column density PDFs 




Fig. 7. Same as Figure [4] but the time- and projection-averaged 
logarithmic column density PDFs of 77 = ln(E, / (2,)) are shown. 
Ej and (E,) denote the column density and the mean column 
density integrated along the i = x, y, z principal axes respec- 
tively. As for the volumetric PDFs of Fig. [4] the standard de- 
viation of the column density PDF obtained from compressive 
forcing is roughly three times larger than from solenoidal forc- 
ing (see Table [3]). The inset in the upper right corner shows the 
PDFs of column density computed in z-projection at a fixed time 
t — 2 T, corresponding to the snapshots shown in Figure [2] The 
density dispersions computed for these instantaneous PDFs are 
crj; = 0.49 and cr n = 0.45 for solenoidal forcing, and cr^ = 1.34 
and o-jj = 1.56 for compressive forcing. Although these distribu- 
tions are quite noisy, the influence of the forcing is still clearly 
discernible. Thus, by studying instantaneous column density 
PDFs, which are accessible to observations, one should be able 
to distinguish solenoidal from compressive forcing. 



Our forcing analysis provides the solution to this apparent dif- 
ference, which l ies at the heart of the d isag reement of the PDF 
data ana lysed in Mac Lo w et all d2005l) andlGlover & Mac Low! 
d2007bl) with the model bvlPassot & Vazquez-Semadenil (119981) . 
iPassot & Vazquez-Semadenil (1 19981) used ID models. In ID, 
only compressive forcing is pos sible, because n o tran sverse 
waves can exist. In contrast, 

|m ac Low et al. I J2001 and 
iGlover & Mac Lowl (l2007bh used a mixture of solenoidal and 
compressive forcing in 3D. In this section, we show that the pa- 
rameter b in both equations (ITST ) and ( fT9b is a function of the 
forcing parameter (. Indeed, using th e relati on cr s = bAi anal- 
ysed in Passot & Vazquez-Semadenil (1998]), but with a lower 
proportionality constant (b = 0.5 in contrast t o b = 1) gives 



a very good representation of the PDF data in iMac Low et al.l 
(2005, Fig. 8). Thus, an investigation of the parameters that con- 
trol b seems necessary and important. 

Moreover, relations ([T8l > and ([T9l are key ingredients 
for the analytical models of the ste llar initial mass func- 
tion by Padoan & N ordlundl d2002l) and iHennebelle & Chabrieil 



d2008l). as well as for t he s tar formation r ate model by 
Krumholz & McKee (2005) and fkrumholz etal. (2009) and for 
the star formation efficiency model by Elm egreenl (120081) . In all 
these models, b is assumed to be 0.5, which is an empirical re- 
sult fr om magnetohydrodynamic a l simulatio ns by Padoan et al.l 
On the other hand, IPassot & Vazquez-Semadenl 
(199 81) found b » 1 f rom ID hydrodynamical simulations. 
Fede rrath et al.l (l2008bh resolved this disagreement betw een 
Padoan et al.ld 19971) and lPassot & Vazquez- Semadenil d 1998b by 



showing that b is a function of the ratio £ e [0, 1] of compres- 
sive to solenoi d al mod es of the turbulence forcing. However, 
iFederrafh et al.l (1200 8b) only tested the two limiting cases of 
purely solenoidal forcing = 1) and purely compressive forc- 
ing (£ = 0). They approximated the regime of mixtures with a 
heuristic model, which had a linear dependence of b on (: 



1 + 



D 



If 



1 - §f, for D = 3 
1 - U, for D - 2 



(20) 



1, 



for D = 1 



Here, we refine this model based on eleven additional sim- 
ulations with f = [0, 1] separated by A£ = 0.1 for RMS Mach 
numbers of 5 in 2D and 3D. These simulations allow us to test 
eleven different mixtures of forcing controlled by the param- 
eter £ (see eq. [9]). The models were run at a numerical res- 
olution of 256 3 and 1024 2 grid points in 3D and 2D, respec- 
tively. We use a lower resolution in 3D, because using our stan- 
dard resolution of 1024 3 would be too computationally intensive. 
However, as shown in § 13.41 the standard deviation of the density 
is fairly well reproduc ed at 256 3 , as is the RMS Mach number M 
dFederrath et al.|[2009l Fig. 2), which allows a reasonably accu- 
rate determination of b. The results are plotted in Figure[8]as dia- 
monds for 3D (top panel) and 2D (bottom panel). We used equa- 
tion dT8l to measure b, because unlike equation ( fT9] >, this version 
of the standard deviation-Mach number relation does not rest on 
the assumption of a log-normal PDF. In fact, if equation ( fl9l i was 
used to derive b for models with £ < 0.5, b would be overesti- 
mated significantly (by up to an order of magnitude for £ = 0), 
because the deviations from the perfect log-normal distributio n 
are stronger for £ < 0.5 (cf. §[321 see also lSchmidt et al.ld2009l) ). 

Figure [8] shows that the dependence of b on £ is non-linear. 
For 3D turbulence the parameter b increases smoothly from 
b « 1/3 for ( = 1 to b * 1 for f = 0, and for 2D tur- 
bulence from b x 1/2 for £ = 1 to b « 1 for £ — 0. 
However, there is an apparent break at £ as 0.5, which repre- 
sents the natural forcing mixture used in many previous stud- 
ies. For f > 0.5 the ^-parameter remains close to the value 
obtained for purely solenoidal forcing, i.e. b * 0.3 - 0.4 in 
3D and b « 0.5 in 2D. The flat part of the data in Figure [8] 
for f > 0.5 explain s why in previous studies with a natural 
forcing mixture (e.g., Mac Low et al. 1998; Kles sen et al.l ]2000; 
iLi et al.ll2003UKritsuk et alj2007l:lGlover et alj2009l) . the turbu- 
lence statistics were close to the purely solenoidal forcing case 
(e.g.JPadoan et alll997l:IStone et alJl998blBoldvrev et al.12 002: 
Pado an & Nordlundl2002l:|Kowal et al.ll2007l;lLemaster & Stond 
2008; Bu rkhart et al.ll2009l) . In contrast, b increases much more 
strongly for £ < 0.5, u ntil it reaches b » 1 for purely 
compressive forcing ( e.g., Passo t~& Vazquez-Semadenil [1998; 
iFederrath et alJl2008btlSchmidt et alj|2009l) . 

Equation ( 120b thus needs to be refined to account for the non- 
linear dependence of b on the forcing. Moreover, equation d20b 
was based on the analytic expression of the forcing parameter 
(T (cf- § 12.11 ). However, the numerical estimate of b depends on 
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Fig. 8. Diamonds: The proportionality parameter b in the den- 
sity dispersion-Mach number relation, eq. ([T8l . computed as 
b - o-pl {{p) M) for eleven 3D models at numerical resolution 
of 256 3 grid cells (top panel) and eleven 2D models at numeri- 
cal resolution of 1024 2 grid cells (bottom panel), ranging from 
purely compressive forcing (( - 0) to purely solenoidal forc- 
ing = 1). The parameter b decreases smoothly from b m 1 
for compressive forcing to b * 1/3 in 3D and b m 1/2 in 2D 
for solenoidal forcing. Stars: Ratio (¥) = E\ ong /E tot of longitu- 
dinal to total power in the velocity power spectrum (see § 16.11 ). 
This quantity provides a measure for the relative amount of com- 
pression induced by the turbulent velocity field, and appears to 
be correlated with the standard deviation of the density PDF. 
Squares: Same as stars, but multiplied by the geometrical fac- 
tor Vd with D — 3 for the three-dimensional case and D — 2 
for the two-dimensional case. The quantity Vd ( v f) provides a 
good numerical estimate of the PDF parameter b. The dashed 
lines show model fits using equation d23l for D — 3 (top panel) 
and D — 2 (bottom panel). 



how well the code can actually induce compression through the 
build-up of divergence in the velocity field. Thus, different codes 
can produce slightly different values of b for the same forcing 
parameter This is because of the varying efficiency of codes 
to convert the energy pr ovided by a given forcing into actual ve- 
locity fluctuations (e.g. lKitsionas et al1l2009h IPrice & Federrathl 
2010). To construct a refined model for b that does not directly 
rest on the analytic forcing parameter f and that accounts for 
the non-linear dependence on the forcing, we recall that b is a 
normalised measure of compression. Compression is caused by 
converging flows and shocks, which have a finite magnitude of 
velocity divergence. A normalised measure of compression is 
thus also provided by dividing the power in longitudinal modes 
of the velocity field by the total power of all modes in the veloc- 
ity field, 



-Elone 



We therefore expect a dependence of b on { y ¥). 

Figure|8]shows (*F) as a function of £ (plotted as stars) for 3D 
and 2D turbulence. It is indeed correlated with b, however, {*¥) 
is less than b by a factor of roughly V3 in 3D and V2 in 2D. The 
squares in Figure[8]show V3 <*P> in 3D and V2 <¥) in 2D, which 
seems to provide a good estimate of b. The factor V3 is a geo- 
metrical factor for 3D turbulence (the diagonal in a cube of size 
unity). It is V2 in 2D turbulence (the diagonal in a square of size 
unity), and VT in ID. The latter in particular is trivial, because 
in ID only longitudinal modes can exist, and thus VI ( V P) = 1 
for any value of £ (cf. Fig. Q]). The larger geometrical factors 
in 2D and 3D account for the fact that the longitudinal velocity 
fluctuations, which induce compression occupy only one of the 
available spatial directions (two in 2D and three in 3D) on av- 
erage. For the general case of supersonic turbulence in D = 1,2 
and 3 dimensions, these ideas lead to 



b = VZ>0F> 



(22) 



which is solely based on the ratio of the power in longitudinal 
modes in the velocity field to the total power of all modes in the 
velocity field, ( V F). 

In addition to the refined model based on the compressive 
ratio {*¥) in equation (l22l . we provide a fit function for b based 
on the forcing parameter (. The dashed lines in Figure [8] show 



b(0 = — + — 
h D D 



1 ( Fi oas (Q 
F lol (0 



(23) 



(21) 



The forcing ratio F\ ong /F lot is given by equation [9] The first 
summand in equation ( 1231 is the expected ratio of longitudi- 
nal modes (compression) in a supersonic turbulent medium for 
a purely solenoidal forcing, i.e. a forcing that does not directly 
induce compression. The second summand is the contribution 
to the compression directly induced by the forcing. The model 
equation ( 1231 is similar to equation ( 1201 . but with a non-linear 
dependence of b on the forcing parameter ( . 

We suggest that t he dependen c e of b on the forcing solves 
a puzzle reported by iPineda et alj d2008l) . They provided mea- 
surements of velocity dispersions and 12 CO excitation tempera- 
tures for the six subregions in the Perseus MC. The molecular 
excitation temperatures serve as a guide for the actual gas tem- 
perature, from which the sound speed can be estimated. From 
these values, the local RMS Mach numbers are computed as the 
ratio of the loc a l velo city dispersion to the lo cal sound speed. 
iGoodman et al.1 d2009h and lPineda et all (120081) pointed out that 
there is clearly no correlation of the form suggested by equa- 
tion ( fT9b for a fixed parameter b across the investigated subre- 
gions in the Perseus MC. For instance, the Shell region exhibits 
an intermediate to small velocity dispersion derived from I2 CO 
and 13 CO observations, while its density dispersion is the largest 
in the Perseus MC. This provides additional support to our sug- 
gestion that the Shell in Perseus is dominated by compressive 
turbulence forcing for which b takes a higher value compared 
to solenoidal forcing. The apparent la ck of density disper sion- 
Mach number correl ation reported by IPineda et al. (2008]) and 
IGoodman et al.l (120091) for a fixed parameter b can thus be ex- 
plained, because b is in fact not fixed across different subregions 
in the Perseus MC. 

We plan to measure b in different regions of the ISM in future 
studies. However, the main problem in a quantitative analysis of 
equation (fT8l with observational data is that the column density 
dispersion is typically smaller than the 3D density dispersion 
(compare Tab. Q] and Tab. [3). The relation between the column 
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density PDF and the volumetric density PDF is non-trivial and 
depends on whether the column density tracer is optically thin 
or optical ly thick and on th e scale of the turbulence driving. 
However, Brunt et al. (2010) developed a promising technique 
to estimate the 3D density variance from 2D observations with 
an accuracy of about 10%. 

4. Intermittency 

Intermittency manifests itself in 

i) non-Gaussian (often exponential) wings of PDFs of quan- 
tities involving density and/or velocity, its derivatives (e.g., 
vorticity) and combinations of density and velocity (e.g., 
p'/ 2 v and p ' 3 v as discussed in AppendixlAl. 

if) anomalous scaling of th e higher-order structu re functions of 
the velocity field (e .g.. lAnselmet etai]ll984l) and centroi d 
velocity increments dLis et al.lll996tlHilv-Blant et al.ll2008l) . 
and 

Hi) coherent structures of in t ense vorticity (V x v ) (see 
I Vincent &"M eneguzz illl99U iMoisv & Jimenezll2004l for re- 
sults of incompressible turbulence), and of strong shocks and 
rarefaction waves (V ■ v). 

Filamentary coherent structures of vorticity (intermittency 
item Hi) are indeed observed in our two supersonic models. 
In Figure [2] (middle panel), we show the projected vorticity 
for solenoidal and compressive forcing, respectively. Most of 
the filaments of high vorticity coincide with the positions of 
shocks and therefore also with high density and negative di- 
vergence in the velocity field (Figure [2] bottom panel). This 
is further more inferred from ob servations of the Ursa Majoris 
Cloud by Falg arone et al.l (I 19941) and is consistent with the re- 
sul ts of weakly com pressibl e decaying t urbule nce experiments 
by iPorter etafl (1 1 992al) and lPorter et all d 19941) . who concluded 
that intense vorticity is typically associated with intermittency. 

4. 1 . The probability distribution of centroid velocity 
increments 

Since there is evidence of filamentary coherent structures in 
the vorticity (intermittency item Hi) of our models, and because 
there is additional evidence of non-Gaussian tails in the density 
PDFs (intermittency item i) discussed in § [3] we now proceed 
to examine the PDFs and the scaling of centroid velocity incre- 
ments (intermittency item ii) to assess the strength of the inter- 
mittency. We compare centroid velocity increments (CVIs) for 
solenoidal and compressive forcing and discuss the interpreta- 
tion of obs ervations based on that comparison. Following the 
analysis by iLis et alJ dl996l), who discuss CVI s computed for 
the turbulence simulation bv lPorter et al.l {T994), and following 
the CVI analysis of the Polaris Flare and of the Taurus MC by 
iHilv-Blant et al.l (I2008I) . the centroid velocity increment is de- 
fined as 



5C c (r) = (C(r) - C(r + 0) 



(24) 



where the angle average ( ) is computed over all possible di- 
rections of the vector t in the plane perpendicular to the line-of- 
sight. Thus, 6Cc(r) only depends on the norm of the lag vector 
I = \t\, which separates two points r = (x,y) and r + ( in the 
plane of the sky (x,y). The normalised centroid velocity, C(r) in 
equation (l24l) is defined as 



C(r) 



j p(r,z)v z (r,z)dz 
J p(r,z)dz 



(25) 



The variable v z (r, z) denotes the line-of-sight velocity in z- 
direction. We have however computed C(r) separately along 
each of the three principal lines-of-sight x, y and z of our 
Cartesian domain in order to examine the effects of varying the 
projection. Also note that we have comput ed normalised cen- 
troid veloci ties (Laz arian & Esquivel 2003), since we want to 
compare to IHilv-Blant et alJ (I2008I) . Another point to mention 
here is that the centroid velocities, C(r) are typically computed 
using an intensity weighting instead of a density weighting. This 
is because the gas density cannot be measured directly, whereas 
the emission intensity is accessible to observations. By using 
density weighting we implicitly assume optically thin emission. 
For opticall y thick emission , uniform weighting would be more 
appropriate dLis et alJ[T9 96). 

Figure|9]shows the PDFs of 5C((r) computed for varying lag 
I in units of the numerical cell si z e A = L/1024. They should 
be compared to IHilv-Blant etal] d2008l Fig. 4-6). The PDFs 
are mainly Gaussian for large lags, whereas for smaller separa- 
tions, they develop exponential tails, indicating intermittent be- 
haviour. Th is result i s consi stent with the numerical simulation 
analysed bv lLis et al.l dl996l) . and with observati ons of the p Op h 
Cloud, the Orion B and the Polaris Flare bv ILis et al.l dl998l) . 
lMiesch etalJ (fl999) and lHilv-Blant et al.l (l2008l h respectively. 

Following the analysis bv lHilv-Blant et ail (120081) . we com- 
puted the kurtosis 7C of the PDFs of CVIs using the definition 
in equations (1131 . Note that 7C = 3 corresponds to a Gaussian 
distribution, and 7C = 6 corresponds to an exponential function. 
The kurtosis of the CVI PDFs is shown in FigurefTOlas a function 
of spa tial lag I, and can be directly compared to Hilv-Blant et al.l 
(2008, Fig. 7). Both forcing types exhibit nearly Gaussian val- 
ues of the kurtosis at lags £ > 100A. On the other hand, for 
( < 100A, both forcing types produce non-Gaussian PDFs. 
Solenoidal forcing approaches the exponential value 'K - 6 for 
I < 10 A. Compressive forcing yields exponential values already 
for lags i « 40A, while solenoidal forcing has K ~ 4 on these 
scales. This indicates stronger intermittency in the case of com- 
pressive forcing. For ( < 30A, compressive forcing yields even 
super-exponential values of 'K. For both solenoidal and com- 
pressive forcings, we show later in § [6] that I < 30A is in the 
dissipation range for numerical turbulence. Compressive veloc- 
ity modes dominate in this regime (see Fig.lT4l). which may re- 
sult artificially in extreme intermittency. For I ss 30A, compres- 
sive forcing gives 7C = 6.0 ± 1.0, which is roughly 35% larger 
than the Polaris Flare observations at their resolution limit. The 
solenoidal case on the other hand gives 7C ~ 4.3 ± 0.5, which 
is in very goo d agreement w i th the IRAM and KOSMA data 
discussed by Hilv- Blant et alJ (120081 Fig. 7). Depending on the 
actual lag used for the comparison, both solenoidal and com- 
pressive forcing seem to be consistent with the observations. 
However, it should be noted that the lags cannot be easily com- 
pared for the real clouds and the simulations, because simulated 
and observed fields have different spatial resolution. Moreover, 
the simulated fields have periodic boundaries, while the true 
fields don't. Nevertheless, the similarity of the observed and the 
numerically simulated CVIs indicates that turbulence intermit- 
tency plays an important role in both our simulations and in real 
molecular clouds. 

The Polaris Flare has a very low star formation rate and is 
therefore appropriate for studying the statistics of interstellar su- 
personic turbulence without contamination by internal energy 
sources. In contrast, the Taurus MC is actively forming stars. 
Against our expectations, the Taurus MC data display very weak 
intermittent beha viour and the kur t osis re mains at the Gaussian 
values 7C ~ 3 in IHilv-Blant et alJ d2008l Fig. 7). However, the 
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Fig. 9. PDFs of centroid velocity increments, computed using equations (1241 and d25T l are shown as a function of lag I in units of grid 
cells A = L/1024 for solenoidal forcing (left) and compressive forcing (right). The PDFs are very close to Gaussi an distributions 
for lo ng lags, whereas for short lags, they develop exponential tails, which is a manifestation of intermittency (e.g.. lHilv-Blant et al.l 
2008, and references therein). 



i± 10- 




l / A 



Fig. 10. Kurtosis 7C of the PDFs of centroid velocity increments 
shown in Fig. [9] as a function of the lag I in units of grid cells 
A = L/1024 for solenoidal and compressive forcing. Note that 
a kurtosis value of 3 (horizontal dot-dashed line) corresponds to 
the value for a Gaussian distribution. Non-Gaussian values of the 
kurtosis are obtained for i < 100A. The error bars contain both 
snapshot-to-snapshot variations as well as the variations between 
centroid velocity increments computed by integration along the 
x, y and z axes. This figure can be compared to obser vations of 
the Po laris Flare and Taurus MC (see Fig. 7 of iHilv-Blant et all 
2008). 



Taurus field studied by IHilv -Blant et al ] (12008b is located far 
from star-forming regions in a translucent part of the Taurus MC 
(E. Falgarone 2009, private communication). This may explain 
why the Taurus field displays only very weak intermittency. It 
would be interesting to repeat the analysis of centroid veloc- 
ity increments for regions of confirmed star formation, includ- 
ing regions with winds, outflows and ionisation feedback from 
young stellar objects to see whether these regions indeed display 
stronger intermittency. 



4.2. The structure function scaling of centroid velocity 
increments 

In this section, we discuss the scaling of the p th order structure 
function of CVIs, defined as 



CVISF p tf) = {\6Ct(r)\ p ) r 



(26) 



We have averaged over a large enough sample of indepen- 
dent increments 8C((r) that increasing the sample size produced 
no change in the value of CVISF p (l) for p < 6, which is demon- 
strated in Appendix 151 FigureQ~T|shows the CVI structure func- 
tions for solenoidal and compressive forcing. The CVI structure 
functions were fit to power laws of the form 



CVISF p (f) oc eh 



(27) 



within the iner t ial ran ge 2 , defined equivalently to the study in 
iFederrafh et al.l (120091) . The value for each power-law exponent 
is indicated in Figure fTTI and summarised in Table|4] 

For a direct comparison of C VI structure functions with the 
study bv lHilv-Blant et all (120081, Fig. 8), we apply the extended 
self-similarity (ESS) hypothesis (IBenzi et al.l 19931) . which states 
that the inertial range scaling may be extended beyond the iner- 
tial range, such that power-law fits can be applied over a larger 
dynamic range. The ESS hypothesis is used by plotti ng the p fh 
order CVISF p (£) against the 3rd order CVISF 3 (f) (IBenzi et al.l 
119931) . These plots are shown in Figure [12] Indeed, the scaling 
range is drastically increased using ESS. All ESS data points are 
consistent with a single power law for each CVI structure func- 
tion order p < 6. We summarise the scaling exponents with and 
without using the ESS hypothesis in Table|4] 

Table |4] furthermore prov ides the ESS scaling ex ponents ob- 
tained for the Polaris Flare dHilv-Blant etal.ll2008l Tab. 3), as 



B y its f ormal definition for incompressible turbulence studies (e.g., 



Frisc the inertial range is the range of scales for which the tur- 

bulence statistics are not directly influenced by the forcing acting on 
scales larger than the inertial range, and not directly influenced by the 
viscosity acting on scales smaller than the inertial range. The inertial 
range is typically very small in numerical experiments, because of the 
high numerical viscosity caused by the discretisation scheme, given the 
resolutions achievable with current computer technology (see also § [6] 
and Appendix let. 
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Fig. 11. Scaling of the structure functions of centroid velocity increments defined in equation d26l i for solenoidal forcing (left) and 
compressive forcing (right) up to the 6th order. Scaling exponents obtained using power-law fits following equation ( f27T > within the 
inertial range are indicated in the figures and summarised in Tab. [4] 

Table 4. Scaling of the structure functions of centroid velocity increments. 
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CVI SFs(1024 3 sol) 
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CVI SFs(1024 3 comp) 
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Relative Scaling Exponents 
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CVI SFs using ESS" (1024 3 


sol) 


0.36 


0.70 


1.00 


1.27 


1.51 


1.72 


CVI SFs using ESS" (1024 3 


comp) 


0.38 


0.72 


1.00 


1.23 


1.41 


1.56 


Polaris Flare* 




0.37 


0.70 


1.00 


1.27 


1.53 


1.77 


Polaris Flare c 




0.38 


0.71 


1.00 


1.28 


1.54 


1.80 


Intermittency Model SL94 rf 




0.36 


0.70 


1.00 


1.28 


1.54 


1.78 


Intermittency Model B02 f 




0.42 


0.74 


1.00 


1.21 


1.40 


1.56 



" Using extended self-similarity (ESS) dBenzi et al.lll993l). 

* Measurement of CVI structure functions by Hilv-Bl ant et alJ d2008h . 

c Measurement of CVI structure functions bv lHilv-Blant et al.1 J2008r) using 12 CO(2-l) data bv lBensch et al.1 d200lh . 

d Intermittency model = p/9 + c(l - (1 - 2/(3C)) p/3 ) defined in eq. {28]! using a fractal co-dimension C = 2 dShe & Levequell 19941) . which 
corresponds to filamentary structures (D = 1). 

e Same as d , but for co-dimension C = 1 (Boldyrev 2002; Boldvrev et al. 2002) corresponding to sheet-like structures (D = 2). 



well as the scaling exponents obtained from intermittency mod- 
els of the structure function scaling exponents 

~ Zp P (I 2 \ p/3 \ 

w °H' -('-3c) I <2S) 

bv IShe & Leveouel (fl99l (C = 2) and iBoldvrevI J200l (C = 
1). In these models, the fractal co-dimension C is related to the 
fractal dime nsion of the mos t inte rmittent structures D by C — 
3-D. The IShe & Leveq ue (1994) model assumes ID vortex 
filam ents as the mo st intermittent structures (T> - 1), whereas 
the IBoldvrevI (j2002) model assumes sheet-like structures with 
D = 2. 

For solenoidal forcing, the scaling of the CVI structure func- 
tions using ESS is very similar to the IShe & Levequel (11994 
model. This model is appropriate for incompressible turbulence, 
for whi ch the most intermitte nt structures are expected to be fila- 
ments (IShe & Levequd[T994l D = 1). Interestingly, their model 
seems t o be consistent with the measurements in the Polaris 
Flare by iHilv-Blant et al.1 (120081) and with our solenoidal forc- 
ing case. In contrast, the scaling exponents derived for compres- 
siv e forcing a re bet ter consistent with the intermittency model 
by IBoldvrevI (120021 D = 2). This direct comparison indicates 



that tu rbulence in the Polaris Flare observed bv lHilv-Blant et ail 
(2008) behaves like solenoidally forced turbulence. However, it 
does not imply that turbulence in the Polaris Flare is close to 
incompressible, since our numerical models are clearly super- 
sonic in the inertial range (see § |SJ. It rather means that CVI 
scaling is different from the absolut e scaling exponents fo llow- 
ing from the int ermittency models bv lShe & Levequel (1 19941) and 
IBoldvrevI <|2002). This is mainly because of two reasons: First, 

these models do not account for density fluctuations (see how- 

1 1 1 1 

ever[Schmidt et al. 2008), and second, CVIs are 2D projections 
of the 3D turbulence. The statistics derived from CVIs is a con- 
volution of density and velocity stat i stics p rojec ted onto a 2D 
plane . As shown bv lOssenkopf et aD (|2006) and lEsquivel et al.1 
d2007h . CVI statistics differ significantly from pure velocity 
statistics, if the ratio of density dispersion to mean density is 
high. This is usually the case in supersonic flows, and is also 
the case for both our numerical experiments (see Table [TJ. 
It explains the difference between the structure functions de- 
rived from the pure velocity statistics compar ed to convolved 
velocity-density statistics (Schmi dt et al.ll20 08). The deviations 
from the Kolmogo rovl (|1941|) s c aling [£ p - p/3) for the 3D 
data analysed in Sch midt et al.l d2008l) are significantly larger 
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Fig. 12. Same as Fig.[TT] but using the extended self-similarity hypothesis (Benzi et al 



the scaling exponents of centroid velocity increments with the study by Hilv-Blant et al 
(seeTab.il 



[1991 allowing for a direct comparison of 



(2008) for the Polaris Flare and Taurus MC 



than those derived via CVI in 2D, revealing a significant loss 
in the si gnatures of inter mittency in the projected CVI data 
(see also iBrunt et~aT]|2003h iBrunt & Mac Lowll2004 for a dis- 
cussion of projection effects). This also means that direct tests 
of the theoretical models will be very difficult to achieve, un- 
less a means of relating the CVI-based moments to the 3D mo- 
ments is developed. Moreover, the fractal dimen sion of struc- 
tures changes in a non- trivial way upon projecti on (IStutzki et alJ 
U998tlSanchez et alj2005tlFed errath et al.l20 09), which severely 
limits the compa rison of CVI st a tistics with the 3D inte rmit- 
tency models by IShe & Levequd (11994 . iBoldvrevi d2002h and 
ISchmidtet al] (120081) . 

Nevertheless, a direct comparison of CVI structure function 
scaling obtained in numerical experiments and observations can 
provide useful information to distinguish between different pa- 
rameters of the turbulence, as for instance different turbulence 
forcings. 



5. Principal component analysis 

Principal compon e nt an alysis (PCA) is a multivariate tool 
(iMurtagh & Heckl Il987l) introduced by iHever & Schloerbl 
(119971) for measuring the scaling of interstellar turbulence. It 
has been used for studying the structure and scaling in sev- 
eral molecular cloud regio n s, simulations and synthetic im- 
ages (IBrunt & Heyeril20"0 2a b; Bru nt et al.ll2003t IHever & Bruntl 
2004: lHever et al.l l2006). PCA can be used to characterise struc- 
ture on different scales. For best comparison with observa- 
tions, we choose to work in position-position-velocity (PPV) 
space. Since our simulation data are typically stored in position- 
position-position (PPP) space, we transformed our PPP cubes 
into PPV space prior to PCA. As for the CVIs discussed in the 
previous section, we use the approximation of optically thin ra- 
diative transfer to derive radiation intensity. This means that we 
essentially assume that the emission is proportional to the gas 
density. The PPV data therefore represent a simulated measured 
intensity T(xi,yi,v z j) = Tjj at spatial position r, = (x,',y,) and 
spectral position v z j. The indices i and j thus represent the spa- 
tial and spectral coordinates respe ctively. A detailed descri ption 
of the PCA technique i s given bv lHever & Schloerbl (Il997l) and 
IBrunt & Heverl d2002al) . The most important steps necessary to 
derive the characteristic length scales and corresponding veloc- 



ity scales using PCA are described below. First, the covariance 
matrix 



1 



N x Ny 



Tij Tik 



(29) 



is constructed by summation over all spatial points N x and N y . 
Solving the eigenvalue equation 

C/«®=A®« (0 (30) 

yields the Ith eigenvalue and the Ith eigenvector it® of the 
covariance matrix. The subsequent projection 



i ik k 



(31) 



,</) 



onto the eigenvectors yields the Ith eigenimage IV 
Autocorrelation functions (ACFs) are then computed for 
each of the eigenimages and eigenvectors. The spatial scale on 
which the two-dimensional ACF of the I th eigenimage falls off 
by 1/e defines the Ith characteristic spatial scale. Following 
the same procedure, the corresponding characteristic velocity 
scale is determined from the ACF of the I th eigenvector, which 
contains the spectral information. 

Figure Q~3] shows our time- and projection-averaged set of 
spatial and velocity scales obtained with PCA. We have fitted 
power laws to the PCA data, which yielded PCA scaling expo- 
nents apcA for solenoidal and compressive forcing respectively. 
For solenoidal forcing we find ckpca = 0.66 ± 0.05 and for com- 
pressive forcing we find «pca = 0.76 ± 0.09 (see Table [5]). The 
different PCA slopes apcA derived for solenoidal and compres- 
sive forcing suggest that by using PCA, differences in the mix- 
ture of transverse and longitudinal modes of the velocity field 
can be detected. However, the difference between solenoidal and 
co mpressive forcing is only at the 1-cr level. 

IHever et all (120061) applied PCA to the Rosette MC and to 
G216-2.5 (Maddalenas's Cloud). These two clouds are quite dif- 
ferent in dynamical and evolutionary state, althoug h they ex- 
hibit roughly the same turbulence Mach number. IHever et aH 
(2006) measured the Mach number Mi pc ~ 4-5 on a scale of 
1 pc for both clouds. The Rosette MC exhibits confirmed mas- 
sive star formation, whereas G216-2.5 has a low star formation 
rate, similar to the Pola ris Flare discussed in the previous sec- 
tion. [He^eretaj] (|2006) measured PCA slopes for both clouds 



Federrath et al.: Turbulence forcing in simulations and observations 



17 




0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 

log„(e/A) log, («/A) 

Fig. 13. Principal component analysis (PCA) for solenoidal (left) and compressive forcin g (right). The PC A slopes obtained for 
solenoidal and compressive forcings are summarised and compared with observations bv lHever et al.l d2006l) in Table [5] The error 
bars contain the contribution from temporal variations and from three different projections along the x, y and z-axes. The data were 
re-sampled from 1024 3 to 25 6 3 grid points prior to PCA. The re-samp ling speeds up the PCA and has virtually no effect on the 
inertial range scaling (see e.g. JPadoan et al.ll2006l: iFederrath et al.ll2009t) . 

Table 5. Comparison of measured PCA scaling slopes. 



1024 3 sol 1024 3 comp Rosette Zone 1 8 Rosette Zone II* G216-2.5 r 

a pca 0.66 ± 0.05 0.76 ± 0.09 

a pca from 12 CO(1-0) 0.79 ± 0.06 0.66 ± 0.06 0.63 ± 0.04 

a PCA from 13 CO(1-0) 0.86 ± 0.09 0.67 ±0.12 0.56 ± 0.02 



8 PCA bv lHever etail i2006h of the interior of an HII region in the Rosette MC. 
b Same a s but exterior of t he HII region. 

c PCA bv lHever et al. (2006) for G216-2.5 (Maddalenas's Cloud). 



and additionally provided the PCA slopes in two distinct sub- 
regions of the Rosette MC. The first subregion is inside the HII 
region (Zone I) surrounding the massive star cluster NGC 2244 3 , 
while the other subregion is outside of this HII region (Zone II). 
The measured PCA slopes obtained from 12 CO and 13 CO obser- 
vations are summarised in Table [5] together with our estimates 
for solenoidal and compressive forcing. The PCA scaling ex- 
ponent for solenoidal forcing is very close to the PCA scaling 
exponents derived from the 12 CO observations in the G216-2.5 
(»pca = 0.63 ± 0.04) and in Zone II of the Rosette MC (a PC A = 
0.66±0.06). In contrast, the PCA slope derived from 12 CO obser- 
vations in Zone I of the Rosette MC (ffpcA = 0.79 ±0.06) is better 
consistent with our compressive forcing case. This indicates that 
Zone I contains more kinetic energy in compressive modes than 
Zone II an d G2 16-2.5. The corresponding 13 CO observations 
reported in iHever etall (120061) yield slightly larger differences 
between the PCA scaling exponents derived for Zone I on the 
one hand, and Zone II and G216-2.5 on the other hand (see also 
Table 0. This supports the idea that Zone I in the Rosette MC, 
and Zone II as well as G216-2.5 contain quite different amounts 
of compressive modes in the velocity field, which may be the 
result of different turbulence forcing mechanisms, similar to the 



3 The formation of the star cluster XA in the Rosette MC was likely 
triggered by the accumulation of material i n the expanding shell sur - 
rounding the OB star cluster NGC 2244 jWang et al.l [2008. 2009). 
This emphasises the importance of expanding HII regions in triggering 
subsequent star formation by compression of gas in expanding shells 
dElmegreen & Ladalll977h . 



differences obtained in purely solenoidal and compressive forc- 
ings. 

6. Fourier spectra 

6.1. Velocity Fourier spectra 

Fourier spectra of the velocity field E( k) are typically used to dis- 
tingu ish between |K olmogorovl d!941l) turbulence, E(k) oc Ar 5 ^ 3 
and iBurgersI (0.948) turbulence, E(k) oc kr 2 . For highly com- 
pressible, isothermal, supersonic, turbulent flow, it has been 
shown that the inert ial range scaling is c lose to Burgers turbu- 
lenc e. For instance , iRritsuk et al.l (I2007I) found E(k) oc AT 1 95 
and ISchmidt et ail (120091) obtained E(k) oc jH- 87 from high- 
resolution numerical simulations. 

The Fourier spectrum of a quantity provides a measure of 
the scale dependence of this quantity. Velocity Fourier spectra 
are thus defined as 

E(k) dk = i JV ■ v* Ank 2 dk , (32) 

where v" deno tes the Fourier transform of the velocity field (e.g., 
lFrischlfl995l) . The total Fourier spectrum can be separated into 
transverse (k ± v) and longitudinal (k || v) parts by apply- 
ing a Helmholtz decomposition. Note that integrating the trans- 
verse energy spectrum yields the kinetic energy in transverse 
(rotational) modes, while integration of the longitudinal energy 
spectrum yields the kinetic energy in longitudinal (compress- 
ible) modes. Furthermore, by integrating the velocity spectrum 
from k\ to ki, one obtains the kinetic energy content on length 
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scales corresponding to the wavenumber interval [&i>&2]- Since 
the mean velocity is zero in our simulations, integration of the 
total velocity Fourier spectrum E(k) over all wavenumbers yields 
the total variance of velocity fluctuations cr 2 v : 



(33) 



The upper bound of the integral is the cutoff wavenumber k c — N 
for a cubic dataset with A^ 3 data points. Thus, k c = 1024 for our 
standard resolution of 1024 3 grid cells. 

In Figure [14] we show the total velocity Fourier spectra 
E(k) as defined in equation d32b together with its decomposition 
into transverse E tmns and longitudinal E\ mg parts for solenoidal 
and compressive forcing respectively. The prominent signature 
of the different forcings on the main driving scale, k — 2 is 
clearly noticeable: Solenoidal forcing excites mostly transverse 
modes, whereas compressive forcing excites mostly longitudi- 
nal modes in the velocity field at k = 2. However, the forcing 
has direct influence only for 1 < k < 3 (see § 12. U . Further 
down the cascade, the turbulent flow develops its own statis- 
tics as a result of non-linear interactions in the inertial range 
5 < k < 15. We emphasise that this scaling range was chosen 
very carefully, since turbulence simulations will only provide a 
smal l inertial range even at resolutions of 1024 3 grid cells (see, 
e.g., iKlein et al]|2007l : iLemaster & Stondl2009|). This is mainl\ 
caused by the bottleneck phenomenon (e.g., iPorter et al.li 1994 



Dobl er et al.1l2003HHaugen & Brandenb urg 2004; Schmid t et al 
2006; Krit suk et al.ll2007l) . which may slightly affect the Fourier 
spectra in the dissipation range. However, the bottleneck phe- 
nomenon had no significant impact on the turbulence statistics 
in our numerical study for wavenumbers k < 40. This is demon- 
strated in Appendix O where we present the resolution depen- 
dence of the Fourier spectra and the dependence on parameters 
of the PPM numerical scheme. We conclude that the statistical 
quantities derived for wavenumbers k < 40 are not significantly 
affected by the numerical scheme or limited resolution applied 
in the present study. 

We apply power-law fits to the inertial range data with 
the resulting slopes indicated in Figure [14] (top panels). 
Both solenoidal and compressive forcing yield slopes con- 
sistent with size-linewidth relations inferred from ob serva- 

119861: 



tions (e.g., 



iMversI [l98j 



1 arson!) 198 lb 
Solomon et al J 1 1 9871: iFalgarone et al 
1994 1; lOssenkopf & Mac Low! |200: 



Perault et al. 
1992b iMiesch & Ballyl 
iPadoan et al.l I2003T 
Hever&Bruntl 12004 IPadoan et all l2006t lOssenkop fetall 
2008b UHeveretal] 12009). and with the r esults of in dependent 
nume r ical simulations (e.g .,|Klessen et al. 2000; Boldvrev et alJ 
2002; IPadoan et al.l 12004: kritsuk et all 120071: ISchmidt et alJ 
2009). Note that size-linewidth relations of the form <t v oc F 



with scaling exponents y = 0.4-0.5 correspond to Fourier spec- 
tra E(k) oc k~P with scaling exponents in the range j3 = 1.8-2.0, 
because y = (J3 - l)/2. However, it must be emphasised 
that the relation between scaling exponents obtained from 
observational maps of centroid velocities (as discussed in § 14.2b 
and 3D velocity fields from simulations is non-trivial, because 
of projection-smoothing and intensity-weighting. Projection- 
smoothing increases the scaling exponents o f the 2D projection 
of a 3D field such that y 7 n - 73d + 1 12 (e.g.JStutzki et all 199 " 
Bru nt & Mac Low! 120041) . However, [Brunt & M ac Low! d200' 
showed that the effect of projection-smoothing is compensated 
statistically (but not identically) by intensity-weighting of 
observed centroid velocity maps. Thus, our measurements of 
velocity scaling seem consistent with observations. 



It is important to note that the transverse parts E lmns {k) fall 
off more steeply than the longitudinal parts E\ ong (k) for both 
forcing types within the inertial range. For solenoidal forcing, 
we find -Etrans(^) 00 £~ L89 and Ei ong (k) oc AT 179 , and for com- 
pressive forcing, E tmns (k) x k~ 2 03 and E\ ong (k) oc k~ l &1 . This re- 
sult indicates that longitudinal modes can survive down to small 
scales, such that comp ression may not be neglec ted anywhere in 
the turbulent cascade. ILemaster & Stond ([2009, Fig. 9, 10) ob- 
tain E tmns (k) K k~ and Ei ong (k) oc £~ L8 for their hydrodynam- 
ical model with solenoidal forcing at a resolution of 1024 3 grid 
points in the Athena code. This is consistent with our findings 
for the scale dependence of the transverse and longitudinal parts 
and shows that the kinetic energy in longitudinal modes must not 
be neglected within the inertial range. 

In order to quantify the relative importance of compression 
over rotation in the turbulent motions, we present plots of the 
ratio 



^long(^) 

E\on g (k) + E tmns (k) E tot (k) 



(34) 



in the bottom panels of Figure [14] Solenoidal forcing yields 
*P « 1 /3 in the inertial range. We emphasise that the ratio *P« 1/3 
was expected from the fact that compression can only occur in 
one of the three available spatial dimensions on average in the 
case of supersonic flow driven by a purely soleno idal force field 
dElmegreen & Scald l2004 iFederrath et ai1l2008bl) . This is the 
fundamental idea on which the heuristic model of the density 
dispersion-Mach number relation given by equation (l20b was 
based. For compressive forcing, we find fs; 1/2 in the inertial 
range as a result of the direct compression induced by compres- 
sive forcing. Thus, solenoidal and compressive forcing produce 
quite similar amounts of compressive modes in the velocity field 
("P « 1 /3 versus ?»1 /2). This means that even fully compres- 
sive forcing can behave very similar to solenoidal forcing in the 
inertial range, as far as pure velocity statistics are concerned. 
However, we show in the next section that the density statistics 
are very different in the inertial range. The same is true for com- 
bined density-velocity statistics (see AppendixlAl. 

We also note here that the rise of *P at k > 40 for both forcing 
types is a numerical effect, which comes from the discretisation 
of the velocity field onto a grid with finite resolution. This shows 
that energy in rotational modes cannot be accounted for accu- 
rately if vortices are smaller than roughly 30 grid cells in each 
direction, whereas longitudinal modes (i.e. shocks) may still be 
well resolved. As a result, the transverse kinetic energy is un- 
derestimated for k > 40 up to the resolution limit k c = 1024. 
However, the plateau of almost constant *P for k < 40 indicates 
that the discretisation had no significant influence on scales with 
wavenumbers k < 40. The effect of underestimating the trans- 
verse kinetic energy due to the discretisation of fluid variables 
is also observed in the ZEUS-3D simulations by Pavlovski et al.l 
(2006, Fig. 2) for wavenumbers k > 10 at numerical resolu- 
tion of 256 3 grid cells. In Appendix ICl we furthermore demon- 
strate that our results for the Fourier spectra are not affected by 
the specific choice of parameters of the numerical scheme for 
wavenumbers k < 40. 



6.2. Logarithmic density Fourier spectra 

In analogy to the velocity Fourier spectra E(k), we define loga- 
rithmic density fluctuation spectra 



S(k)&k. 



(s-(s))(s-(s)yAnk 2 dk. 



(35) 
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Fig. 14. Top panels: Total, transverse (rotational) and longitudinal (compressible) velocity Fourier spectra E(k) defined in equa- 
tion ( T32T ) and compensated by k 2 for solenoidal (left) and compressive forcing (right). Error bars indicate temporal variations, which 
account for an uncertainty of roughly ±0.05 of all scaling slopes reported for the inertial range 5 < k < 15. The inferred inertial 
range scaling exponents for both solenoidal and compressive forcing are consistent with independent numerical simulations and with 
observations of the size-linewidth relation (see text). Note that the transverse part, E tmns falls off more steeply than the longitudinal 
part, E\ ong for both forcing types in the inertial range. Bottom panels: Ratio of the energy in longitudinal velocity modes E\ ong to the 
total energy in velocity modes E tot = E trdns +E\ ong . For solenoidal forcing, we obtain E\ ong /E tot as 1/3 in th e inertial range (horizontal 
dash-dotted line), beca use compression can only occur in one of the three spatial dimensions on average (Elmegreen & Scalo 2004; 
Fede rrath et al.ll2008bl) . For compressive forcing, this ratio is roughly 1/2, which corresponds to an equipartition of longitudinal 
and transverse velocity modes. Note however that compressive forcing can compress the gas i n all three spatial dime nsions directly, 
whereas solenoidal forcing can only induce compression indirectly through the velocity field dFederrath et a l. 2008b). The excess of 
longitudinal modes at high wavenumbers k > 40 stems from numerical dissipation, which is more effectively dissipating transverse 
than longitudinal modes on small scales due to the discretisation onto a grid. This suggests that roughly 30 grid cells are needed 
to accurately resolve a vorte x, while a shock is typically resolved with roughly 3 grid cells using the piecewise parabolic method 
(IColella & Woodwardl[T9 84). However, for a numerical resolution of 1024 3 grid cells, we find that wavenumbers k < 40 are almost 
unaffected by the discretisation and by the parameters of the numerical scheme (see AppendixO. 



We subtract the mean logarithmic density prior to the Fourier 
transformation such that S (k) is a measure of density fluctuations 
as a function of scale. Therefore, integrating S (k) over all scales 
yields the square of the logarithmic density dispersion cr s 



£ S(k)dk = o- 2 s . (36) 

Furthermore, integrating S(k) over the wavenumber range 
[^1,^2] yields the typical density fluctuations on length scales 
corresponding to this range of scales. 

Figure [15] (left) shows the logarithmic density fluctuation 
spectra S (k) together with the total velocity Fourier spectra E(k) 
in one plot. In contrast to the scaling of the velocity E(k), the 
scaling of S(k) cc k~P is significantly different for solenoidal 
(J3 = 1.56 + 0.05) and compressive forcing (fi = 2.32 ± 0.09) 
in the inertial range. 



7. A-variance of the velocity and density 

The A-variance technique provides a complementary method 
for measuring the scaling exponent of Fourier s pectra in the 
physi cal domain using a wavelet transformation ( Stutzki et al. 
1998). We apply the A-variance t o our simulation data usi ng the 
tool developed and provided by lOssenkopf et all d2008a). This 
tool implements an i mproved version o f the original A-variance 
(IStutzki et alJl!998t iBensch et ai1l200ll) . The A-variance mea- 
sures the amount of structure on a given length scale t by filter- 
ing the dataset q(r) with an up-down-function Q( (typically a 
French-hat or Mexican-hat filter) of size £, and computing the 
variance of the filtered dataset. The A-variance is defined as 

rt(V = \(q(r)*Q(r)Jj , (37) 

where the average is computed over all data points at po- 
sitions r = (x,y,z). The operator * stands for the convo- 
lution. We use the original French-hat filter with a diame- 
ter ratio of 3.0 as in previous studies using the A-variance 
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Fig. 15. Left panel: Fourier spectra of the velocity, E(k) defined in eq. ( 1321 (crosses and diamonds) and Fourier spectra of the 
logarithmic density fluctuations, S (k) defined in eq. (f35T > (triangles and squares) for solenoidal and compressive forcing, respectively. 
Both E(k) and S (k) are compensated by k 2 allowing for a better determination of the inertial range scaling. The density fluctuation 
power spectra differ significantly in the inertial range 5 < k < 15 with S {k) oc k~ l 56 for solenoidal and S (k) oc k~ 232 for compressive 
forcing. The scale on which the density fluctuation spectra from solenoidal and compressive forcing cross each other and where the 
slope obtained in compressive forcing breaks and approaches the shallower slope of the solenoidal forcing case roughly coincides 
with the sonic wavenumber k s (vertical dashed lines) defined in eq. (1381 1. Right panel: Same as left panel, but instead of using Fourier 
spectra to determine the inertial range scaling, we use the A-variance method to derive the scaling slopes in physical space. Note 
that the scaling slop es a obtained with the A-variance technique are related to the slopes f3 of the Fourier spectra by (3 — a + 1 
dStutzkietall l998). Error bars denote \-cr temporal fluctuations. 



techn i que (e.g.. IStutzki et all 1998t iMac Low & O ssenkopf 
2000; Ossenk opf et all l200lt bssenkopf & Mac Lowl l2002t 
Ossenkopf etal.l2006h . 



Figure [15] (right panel) shows that the inertial range scaling 
obtained with the A-variance technique is in very good agree- 
ment with the scaling measured in the Fourier spectra. Note 
that the scaling exponents /3 of Fourier spectra are ideally re- 
lated to the scaling exponents a of the A-variance by a — /3 - 1 
(IStutzki et al.l[T998h . The small deviations from this analytical 
relation are caused by the finite size of the dataset, the re- 
sampling procedure prior to the A-vari ance analysis applied her e 
and the choice of the filter function (Ossenkopf et al. 2008a). 
However, these deviations are on the order of 4% and therefore 
smaller than the average snapshot-to-snapshot variations. 

For the A-variance of the velocity field, cr 2 h (v, I) oc we 
find scaling exponents a = 0.83 + 0.05 for solenoidal forc- 
ing and a = 0.96 ± 0.05 for compressive forcing. This trans- 
lates into size-linewidth relations cr A (v, £) oc p with scaling 
exponents y - a/2. Thus, we find y = 0.42 + 0.03 for 
solen oidal forcing and y - 0.48 ± 0.03 for compressive forc- 
ing. Ossenkopf & Mac Low (2002) found a common power-law 
slope y = 0.5 + 0.04 for the Polaris Flare, ranging over three 
orders of magnitude in length scale from about 50 pc down to 
roughly 0.05 pc. This scaling exponent is roughly consistent with 
both our solenoidal and compressive forcing data, but slightly 
better consistent with co mpressive forcing. Note that t he cen- 
troid velocity analysis by Ossenkopf & Mac Low! (120021) is also 
subject to the combined effects of projection- smoothing and 
intensity-weighting discussed in iBrunt & Mac Lowl d2004l) and 
discussed in § 16.11 Thus, the comparison of 3D scaling of the ve- 
locity with 2D observations should always be made with the cau- 
tion that projection-smoothing and intensity -weighting roughly 
cance l each other out in a statistical sense dBrunt & Mac Lowl 
l200l . 

We are not aware of any observational study considering the 
scaling of logarithmic intensity. The use of logarithmic density 



is useful in isothermal simulations, because the equations of hy- 
drodynamics, equations (f2]i and (O, are invariant under transfor- 
mations in s — ln(p/(p)). In observations however, the inten- 
sity, T is measured instead of the density, but the intensity can 
be transformed into s' = ln(T/(T)), which gives a normalised 
quantity similar to s = ln(p/(p)). This enables a straightforward 
comparison of simulation and observational data (yet with the 
limitations listed in §|9]l. It is also interesting to look at logarith- 
mic density and intensity scaling, because this scaling parameter 
is used in analytic models of the mass distri bution of cores and 
stars bv lHennebelle & Chabrierl(l2008ll2009h . 

Unlike a logarithmic scaling analysis, the scaling of the 
linear integrated in tens ity, cr^jp, €) oc t? was analysed by 
IStutzki et al] (1 19981) and IBensch et all d200ll) . They found y * 
0.5-0.9 for the Polaris Flare, in good agreement with the scal- 
ing ex ponent y = 0. 8 + 0.1 obtained from solenoidal forc- 
ing in iFederrath et all <|2009). In contrast, the scaling expo- 
nent obtained fo r compressive forcin g is significantly higher 
(y = 1.4 ± 0.3). IBensch' et al.l d2001l) measured scaling expo- 
nents y as 1.0-1.5 in small-scale maps ({ < 0.1 pc) of the Polaris 
Flare and in Perseus/NGC1333, which are consisten t with our 
estimates for compressive forcing ( Fed errath et al.l20"09l, Tab. 1). 
Since both solenoidal and compressive forcings display strong 
intermittency at short lags (see Fig. ITOt. intermittency appears to 
be primarily measurable on scales smaller th an the turbulence in- 
jection scale. Taking together the results by Bens ch et al.l ([2001 ) 
with ours for solenoidal and compressive forcing indicates that 
interstellar turbulence is driven primarily on large scales, poten- 
tially with a significant amo unt of compressiv e modes present 
on the forcing scale (see also lBrun t et alj [2009l) . 

8. The sonic scale 

The velocity Fourier spectra E(k) discussed in § !6.1l can be de- 
scribed as power laws E(k) oc with negative power-law ex- 
ponents, /? > 1 . This means that the typical velocity fluctuations 
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are decreasing when going to smaller scales. The value of the 
integral J k ' E(k) dk over a finite range of wavenumbers with k\ 
as the lower bound and the cutoff wavenumber k c as the upper 
bound therefore becomes smaller with increasing k\. Thus, the 
turbulent flow is expected to change from a supersonic to a sub- 
sonic flow on a certain length scale. This scale separates the su- 
personic regime on large scales, where the velocity fluctuations 
are supersonic from the subsonic regime, which is located on 
smaller scales, where the typical velocity fluctuations are small 
compared to the thermal motions of th e gas. This transition scale 
is called the sonic scale A s . Following ISchmidt et al.l (120091) . the 
corresponding sonic wavenumber k s in Fourier space is defined 
by solving the equation 

£E(k)dk^ l -cl 



(38) 



implicitly for k s . The sonic scale is thus defined as the scale on 
which the mean square velocity fluctuations become comparable 
to the mean square of the sound speed. 

We solved equation (1381 for the sonic wavenumbers k s for 
both the solenoidal and compressive forcing cases. The sonic 
wavenumbers for solenoidal and compressive forcings are in- 
dicted in Figure[l5](left) as vertical dashed lines. We find k s = 26 
for solenoidal forcing and k s = 27 for compressive forcing. The 
corresponding sonic scales A s are also indicated in Figure Q3] 
(right) as vertical dashed lines. 

The Fourier spectra S (k) shown in Figure Q3] (left) and the 
corresponding A-variance curves shown in Figure Q3] (right) for 
solenoidal and compressive forcing cross each other roughly at 
the sonic wavenumber and on the sonic scale, respectively. For 
compressive forcing S (k) is significantly steeper on scales larger 
than the sonic scale (k < k s ) compared to scales k > k s . S(k) 
for compressive forcing approaches the shallower slope of S {k) 
for solenoidal forcing at k « k s . For k > k s there are neither 
significant differences between the density spectra S (k) nor the 
velocity spectra E(k) for solenoidal and compressive forcings. 

The strong break in the logarithmic density fluctuation spec- 
tra S (k) for compressive forcing around k s appears to be linked 
to the transition from supersonic motions on large scales to sub- 
sonic motions on scales smaller than the sonic scale. In order 
to quantify this, we estimated the typical density fluctuations on 

supersonic scales (k < k s ) by evaluating o- 2 (k<k s ) — s S(k) dk. 
We obtain cr s (k < k s ) ~ 1.22 for solenoidal and cr s (k < k s ) ~ 3.05 
for compressive forcing, which is on the order of the logarith- 
mic density dispersions cr s found from the density PDFs (see 
Table Q]). This means that most of the power in density fluctua- 
tions is located on scales larger than the sonic scale. In contrast, 
on scales smaller than the sonic scale the typical density fluc- 
tuations can be estimated by solving cr 2 s (k > k s ) = J k L S(k)dk. 
We obtain cr s (k > k s ) * 0.45 for both types of forcing. This 
shows that density fluctuations on scales below the sonic scale 
are small compared to the typical dens ity fluctuations in the su- 
perso nic regime at k < k s (see also IVazquez-Semadeni et alJ 
120031) . Moreover, Figure Q3] shows that the typical logarithmic 
density fluctuations are similar for both solenoidal and compres- 
sive forcings on scales smaller than the sonic scale. Note that the 
sum of logarithmic density fluctuations on all scales is 



[o J2 s (k<k s ) + o J2 s (k>k s )] -1.30 
for solenoidal forcing and 



[ojikKks) + o- 2 ,(k>k,)f /2 * 3.08 



(39) 



(40) 



for compressive forcing. As expected from equation ( 1361 ), these 
values are in excellent agreement with the total logarithmic den- 
sity dispersions cr s , obtained from the density PDFs shown in 
Tabled] 

A spatial representation of the structures exhibiting sub- 
sonic velocity dispersions is shown in Figure Q~6] (bottom panel). 
These structures are identified in slices through the local Mach 
number M as regions with M < 1 . Figure [16] (top panel) dis- 
plays the corresponding density slices. The density-Mach num- 
ber correlations are quite weak, as expected for isothermal turbu- 
lence (cf. § 13.61 ). However, Figure [5] shows that high-density re- 
gions exhibit lower Mach numbers on average. In real molecular 
clouds, the sonic scale is expected to be l ocated on length scales 
A s » 0.1 pc within factors of a few (e. g., Falgarone et al. 1992; 



Barr anco & Goodmanl 1 998HGoodman et alJl998HSchnee et all 
120071) . Forinstance. lHever et alJd2006l) found ^~0.3-0.4 pc for 
the Rosette MC and A s « 0.1 -0.2 pc for G2 16-2.5. Furthermore, 
the sonic scale may be associated with the transition to coher- 
ent cores ( Goodman et al.|[T998l:lBallesteros-Paredes et a 



Klessen et alJ 2005). Recent simulations of turbulent core forma- 
tion by ISmith et all (120091) also suggest that star-forming cores 
typically exhibit transonic to subsonic velocity dispersions. This 
can be understood if cores form close to the sonic scale in a 
globally supersonic turbulent medium. Figure [16] suggests that 
regions with subsonic velocity dispersions have different shapes 
and sizes for both solenoidal and compressive forcings. The 
movie (available in the A&A online version) shows that these 
structures are transient objec ts, forming and dissolving in the 
turbulent flow (e.g., see also lVazquez-Semadeni et aLll2005l) . If 
we had included self-gravity in the present study, some of these 
regions would have likely collapsed gravitationally, because tur- 
bulent sup port becomes insufficient in some of these subsonic 



bulent support becomes insurricient in 
cores (e.g jMac Low & Klessenl l2004). 



9. Limitations 

As a result of the simplicity of the hydrodynamic simulations 
presented in this paper, comparisons with observational data are 
limited and should be considered with caution. These limitations 
are listed below: 

- We assume an isothermal equation of state, so our mod- 
els are strictly speaking only applicable to molecular gas 
of low enough density to be optically thin to dust cool- 
ing. Variations in the equatio n of state can lead to changes 
in the density statistics (e.g.. [ Passot & Vazauez-Semadenil 
ll998tlLTeTal]l2003tlAudit & Hennebelldl2009h . The results 
of the present study apply primarily to the dense interstellar 
molecular gas for which an isothermal equation of state is an 
adequate approximat i on dWolfire et all 1995b lFerrierdl200 it 
Pavl ovski et al.ll2006t iGlover et alj|2009t) . 

- The numerical resolution of our simulations is lim- 
ited. As shown in Figure [6] the high-density tails 
of the PDFs systematically shift to higher densities 



[y shitt 

(see also [H ennebelle & Aud it] 12001 [Kitsionas et al.l 120091: 
IGlover et alJl2009t iPrice & Federrath! l2010h . However, the 



mean and the dispersions are well converged at the numerical 
resolutions of 256 3 , 512 3 and 1024 3 grid points used in this 
study. The inertial scaling range is very small even at resolu- 
tions of 1024 3 grid cells. However, the systematic difference 
in the inertial range scaling between resolutions of 512 3 and 
1024 3 grid points is less than 3% (see Appendix ICb. which 
is less than the typical temporal variations between different 
realisations of the turbulent velocity and density fields. 
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Fig. 16. z-slices through the local density (top panels) and Mach number fields (bottom panels) at z — and t — IT for solenoidal 
forcing (left), and compressive forcing (right). Regions with subsonic velocity dispersions (Mach < 1) are distinguished from 
regions with supersonic velocity dispersions (Mach> 1) in the colour scheme. The correlation between density and Mach number 
is quite weak. However, as shown in Fig. [5] high-density regions exhibit lower Mach numbers on average. Thus, dense cores might 
naturally exhibit transonic to subsonic velocity dispersions, because their sizes are expect ed to be comparable to the sonic scale. 
The sonic scale may be the transition scale to coherent cores (e.g.. iGoodman et"ai] [T998). Although many of these 'cores' here 
are transient, some of them are dense enough to become gravitationally bound, and accumulate enough mass to decouple from the 
overall supersonic turbulent flow. See the A&A online version for a movie, showing the time evolution of this figure. 



Our simulations adopt periodic boundary conditions. This 
implies that our simulations can only be representative of 
a subpart of a molecular cloud, for which we study turbu- 
lence statistics with high-resolution numerical experiments. 
However, we cannot take account of the boundary effects 
in real molecul ar clouds. Simulations of large-scale collid- 
ing flows (e.g., Heitsch et al. 2006t Vazauez-Semadeni et alJ 
2006; Henneb elle et alJl2008h l Baneriee et al.l l2009) are more 



suitable for studying the boundary effects during the forma- 
tion of molecular clouds. 

We only analysed driven turbulence. However, there is on- 
going debate about whe t her turbulence is driven or decaying 
(e.g., [Stone et alJll998t lMac Lowlll999t iLemaster & Stond 
2008; Offner et al.l2008l) . We are aware of the possibility that 
turbulence may in fact be excited on scales larger than the 



size of molecular clouds (e.g jBrunt et alJl2009h . but may be 
globally decaying (if not replenished by a mechanism act- 
ing on galactic scales). As discussed in § 12.11 this large-scale 
decay can however act as an effective turbulence forcing on 
smaller scales, because kinetic energy is transported from 
large to small scales through the turbulence cascade. 
- Centroid velocity and principal component analysis were 
applied to PPV cubes constructed from the simulated ve- 
locity and density fields assuming optically thin radia- 
tion transfer to estimate the intensity of emission lines. 
This approximation will of course not hold for optically 
thick tracers. A full radiative transfer calculation taking 
account of the level population (e.g ., Keto et alJ 1200 "' 



Steinacker et al. 2006; Pi nte et al.l2009tlHauschildt & Barou 
2009L iBaron et all 12009) of self-consistently formed and 
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evolve d chemical tracer mo lecules (e.g. JGlover & Mac Low! 
12007 allbt iGlover et a i1 l2009h would be needed to advance on 
this issue. 

We neglected magnetic fields. In ord er to test the role of 
magnetic fields in s tar formation (e.g.. lCrutcher et al]|2009t 
iLunttilaet ail [2008). we would have to include the effects 
of magnetic fi elds and ambipolar diffusion . For instance, the 
IMF model bv lPadoan & N ordlund ( 2002) requires magnetic 
fields to explain the present-day mass function, while it is 
still not clear whether magnetic fields are dyna mically im- 
portan t for typical molecular clouds. However, iHever et all 
(2008) showed that magnetohydrodynamic turbulence in the 
Taurus MC may lead to an alignment of flows along the field 
lines. 

The present study did not include the effects of self- 
gravity, because we specifically focus on the pure turbu- 
lence statistics obtained in solenoidal and compressive forc- 
ings. In a follow-up study, we will include self-gravity and 
sink particles . (e.g.. iBate et al.lll995l: iKrumholz et al.ll2004t 



lJappsen etal1l2005UFederrath etal.Tl2010h to study the in- 
fluence of the different forcings on the mass distributions 
of sink particles. First results indicate that the sink parti- 
cle formation rate is at least one order of magnitude higher 
for compressive forcing com pared to solenoidal forcing. 
Vazauez- Semadeni et al.1 (120031) argue that the star forma- 
tion efficiency is mainly controlled by the RMS Mach num- 
ber and the sonic scale of the turbulence (cf. §[8}. However, 
our preliminary results of simulations including self-gravity 
show that the star formation efficiency measured at a given 
time (i.e., the star formation rate) is much higher for com- 
pressive forcing than for solenoidal forcing with the same 
RMS Mach number and sonic scale. This provides additional 
support to our main conclusion that the type of forcing must 
be taken into account in any theory of turbulence-regulated 
star formation. This needs to be investigated in future, high- 
resolution numerical experiments including self-gravity and 
sink particles. 



10. Summary and conclusions 

We presented high-resolution hydrodynamical simulations of 
driven isothermal supersonic turbulence, which showed that the 
structural characteristics of turbulence forcing significantly af- 
fect the density and velocity statistics of turbulent gas (see also 
ISchmidt et al.l i2009). We compared solenoidal (divergence-free) 
forcing with compressive (curl-free) turbulence forcing. Five 
different analysis techniques were used to compare our simu- 
lation data with existing observational data reported in the liter- 
ature: probability density functions (PDFs), centroid velocity in- 
crements, principal component analysis, Fourier spectrum func- 
tions, and A-variances. We find that different regions in the tur- 
bulent ISM exhibit turbulence statistics consistent with different 
combinations of solenoidal and compressive forcing. Varying 
the forcing parameter f e [0, 1] in equation (O, we showed 
that a continuum of turbulence statistics exists between the two 
limiting cases of purely solenoidal = 1) and purely com- 
pressive forcing (£ = 0). For £ > 0.5, turbulence behaves al- 
most like in the case of purely solenoidal forcing, while for 
£ < 0.5, turbulence is highly sensitive to changes in f (cf. [8]). 
Note that £ = 0.5 represents the natural forcing mixture used 
in many previous turbulence simulations. Because the behaviour 
of all forcing mixtures with £ > 0.5 is similar to that of purely 
solenoidal turbulence with £ = 1 (see [SJ, turbulence statistics 
is biased towards finding solenoidal-like values. However, ob- 



servations of regions around massive stars that drive swept-up 
shells into the surrounding medium (e.g., the shell in the Perseus 
MC and in the Rosette MC) seem better consistent with models 
of mainly compressive forcing (£ < 0.5). Note that expanding 
HII regions around massive stars, and supernova explosions typ- 
ically create such swept-up shells, which are considered to be 
important drivers of interstellar t urbulence dMac Low & Klessenl 
l2004tlBreitschwerdt et al.l2009t) 4 . A detailed list of our results is 
provided below: 

1. The standard deviation (dispersion) of the probability dis- 
tribution function (PDF) of the gas density is roughly three 
times larger for compressive forcing than for solenoidal 
forcing. This holds for both the 3D density distribu- 
tions (Figure [4] and Table Q]) and the 2D column den- 
sity distributions (Figure [7] and Table [3]). We extended the 
density dispersion-Mach number relations, equation ( fT8l 
and ( PT9l ) originally investigat e d by Padoan e t al. ( 1997|) and 
iPassot & V azquez- Semadenil (119981) . Based on the varying 
degree of compression obtained by solenoidal and compres- 
sive forcing, we developed a heuristic model for the pro- 
portionally constant b in the density dispersion-Mach num- 
ber re lation, which takes ac count of the forcing parame- 
ter £ dFederrath et alJ l2008bi) . In the case of compressive 
forcing the proportionality co nstant b is close to b » 1, 
which confirms the result by IPassot & Vazquez-Semadenil 
(119981) . In contrast, solenoidal forcing yields b » 1/3, which 
is in excellent agreement with recent independent high- 
resol ution numerical s imulations using solenoidal forcing 
fe.g JBeetz e"taD2008l) . 

2. A parameter study of eleven models with varying forcing 
parameter ( - [0,1], separated by A£ = 0.1 showed that 
the heuristic model given by equation d20b can only serve 
as a first-order approximation to the forcing dependence of 
b (cf. Fig. [8). We showed that b scales with the normalised 
power of compressible modes in the velocity field, (*P). A 
good approximation for b is given by b « Vd (*¥), where 
D = 3 in 3D turbulence. 

3. We compared the density PDFs i n our models with ob - 
servations in the Perseus MC by [Goodman etjri] (120091) . 
iGoodman et al.1 (l2009h obtained the largest density disper- 
sion in all of the Perseus MC within a region that they call 
the Shell region. This Shell surrounds the massive star HD 
278942 suggesting that the Shell is an expanding HII re- 
gion. Swept-up shells represent geometries that can be asso- 
ciated with compressive turbulence forcing, because an ex- 
panding spherically symmetric shell is driven by a fully di- 
vergent velocity field. This may explain why the Shell region 
in the Perseus MC exhibits the largest density dispersion 
a mong all of the subregio ns in the Perseus MC investigated 
by IGoodman et al] |2009). We emphasise that the Shell re- 
gion does not exhibit the highest RMS Mach number, but has 
an intermediate value among the e xamined subregions in the 
Perseu s MC dPineda et al. 2008). Furthermore, as pointed 
out by IGoodman et alJ ( 20091) the density dispersion-Mach 
number relation of the form given by equation (1191 1 for a 
fixed parameter b is not observed for the Perseus MC. This 
apparent contradiction with equation dT9b for a fixed param- 
eter b is resolved, if different turbulence forcing mechanisms 
operate in different subregions of the Perseus MC, such that 
b is a function of the mixture of solenoidal and compressive 
modes f as shown in Figure [8] 



See also Tamburro et al. ( 200^) for an observational study. 
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The turbulent density PDF is a key ingredient for the analyt- 
ical models of the core mass fu nction (CMF) and the stellar 
initial mass function (IMF) by Padoa n & Nordlundl d2002l) 
and iHennebelle & Chabrierl (120081120091), as well as for th e 
star formation rate mo dels by iRrumholz & McKed (f2 005). 
iKrumholz et"ai1 (120091) and lPadoan & Nordlundl (120091). an d 
the star formation efficiency model by lElmegreenl ([2008). 
We showed that the dispersion of the density probability dis- 
tribution is not only a function of the RMS Mach number, 
but also depends on the nature of the turbulence forcing. All 
the analytical models above rely on integrals over the den- 
sity PDF. Since the dispersion of the density PDF is highly 
sensitive to the turbulence forcing, we conclude that star 
formation properties derived in those analytical models are 
strongly affected by the assumed turbulence forcing mecha- 
nism. 

The PDFs p s (s) of the logarithm of the density s = ln(p/ (p)) 
are roughly consistent with log-normal distributions for both 
solenoidal and compressive forcings. However, the distri- 
butions clearly exhibit non-Gaussian higher-order moments, 
which are associated with intermittency. Including higher- 
order corrections represented by skewness and kurtosis is 
absolutely necessary to obtain a good analytic approxima- 
tion for the PDF data, because the constraints of mass con- 
servation (eq. [TlT l and normalisation (eq. \12\ of the PDF 
must always be fulfilled. Even stronger deviations from 
perfect log-normal distributi ons are expected if the gas 
is non-isotherma l (e.g., Passot & Vazquez-Semadeni[ 19981: 
Scal o et al. 1 11998b iLi et al.l l2003l)~m agnetised (e.g.. iLi et 



i et all 
2004; 



2008) or self-gravitating (e.g..lKle ssen 2000 yLi et al 
Federrath et al.ll2008al : lKainulainen et al.ll2009h . which often 
leads to exponential wings or to power-law tails in the PDFs. 
Non-Gaussian wings of the density PDFs are a signature of 
intermittent fluctuations, which we further investigated us- 
ing centroid velocity increments (CVIs). We find strong non- 
Gaussian signatures for small spatial lags £ in the PDFs of 
the CVIs (Figure[9]). These PDFs exhibit values of the kurto- 
sis significantly in excess of that expected for a Gaussian (see 
FigurefTOt . Figure[10]can be compared with lHilv-Blant et al.l 
(2008, Fig. 7), who analysed CVIs in the Taurus MC and 
in the Polaris Flare. The values of the kurtosis 7C measured 
in the Polaris Flare are consistent with exponential values 
(7C = 6) for short spatial lags, which is also compatible 
with the results of solenoidal forcing. In contrast, compres- 
sive forcing yields values of the kurtosis twice as large at 
small lags, which indicates that compressive forcing exhibits 
stronger intermittency. The scaling of the C VI structure func- 
tions supports the conclusion that compressive forcing ex- 
hibits stronger intermittency compared to solenoidal forcing 
(see Figure [12] and Table |4|. The scaling exponents of the 
C VI structure functions obtained for solenoidal forcing are in 
good agreement with the results by iHilv- Blant et al.l OOOlh 
obtained in the Polaris Flare for the CVI structure functions 
up to the 6th order using the extended self-similarity hypoth- 
esis. 

We applied principal component analysis (PCA) to our mod- 
els. A comparison of the PCA scaling exponents ckpca with 
the PCA study in the Rosette MC and in G216-2.5 by 
iHever etal] d2006l) showed that solenoidal forcing is consis- 
tent with the PCA scaling measured in G216-2.5 and with the 
PCA scaling measured in the outside of the HII region (Zone 
II) surrounding the OB star cluster NGC 2244 in the Rosette 
MC. On the other hand, the PCA scaling inside this HII re- 
gion (Zone I) is in good agreement with the PCA scaling ob- 



tained for compressive forcing (Table[5]). Similar to the Shell 
region in the Perseus MC, the HII region in the Rosette MC 
(Zone I) displays signatures of m ainly compressive forcing . 
Recent numerical simulations by Grits chneder et al] (l2009h 
also show that ionisation fronts driven by massive stars can 
efficiently excite compressible modes in the velocity field. 

8. The Fourier spectra of the velocity fluctuations showed that 
they follow power laws in the inertial range with E(k) oc 
r i.86±o.os for so ienoidal forcing and E(k) oc r L94±005 for 
compressive forcing. Both types of forcing are therefore 
compatible with the scaling of velocity fluctuations inferred 
from observations and independent numerical simulations. 
The Fourier spectra of the logarithmic density fluctuations 
scale as S(k) oc fc-i-56±o.05 f or solenoidal forcing and S (k) oc 
£-2.32±o.09 f or compressive forcing in the inertial range. 

9. The inertial range scaling of the velocity and logarithmic 
density fluctuations inferred from the Fourier spectra was 
confirmed using the A-variance technique. 

10. We computed the sonic scale by integrating the velocity 
Fourier spectra. The sonic scale separates supersonic tur- 
bulent fluctuations on large scales from subsonic turbu- 
lent fluctuations on scales smaller than the sonic scale. We 
found a break in the density fluctuation spectrum S {k) for 
compressive forcing roughly located on the sonic scale. 
The typical density fluctuations computed by integration 
of S (k) over scales larger than the sonic scale are consis- 
tent with the logarithmic density dispersions derived from 
the probability density functions for solenoidal and com- 
pressive forcings. On the other hand, the typical density 
fluctuations on scales smaller than the sonic scale are sig- 
nificantly smaller for both forcing types, which may re- 
flect the transition to coherent cores (e.g., iGoodman et al.l 
fl998h . Indeed, observations show that cores typically have 
transonic to subsonic internal velocity dispersions (e.g. , 
Benson & Mversl [19891: [Andre et all 120071; iKirk et all 200" 



Ward-Thompson et al. 2007; Lada et al. 2008; Foster et al 
l2009t iFriesen et al.l l2009t iBeufher & Henningll2009l) . This 
can be understood if cores form near the sonic scale at the 
stagnation points of shocks in a globally supersonic turbu- 
lent ISM (cf. §©. 
11. We found that the correlations between the local densi- 
ties and the local Mach numbers are typically quite weak 
(Figures [5] and [TBI ). However, this weak correlation shows 
that the local Mach number M decreases with increasing 
density as M(p) oc p-° 06 for solenoidal forcing and M(p) oc 
p-o.05 f or compressive forcing for densities above the mean 
density. This means that dense gas tends to have smaller ve- 
locity dispersions on average, consistent with observations 
of dense protostellar cores. 
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Appendix A: Fourier spectra and A-variance scaling 
of the combined quantities p 1/2 v and p 1/3 v 

In this section we present the Fourier spectra and A-variance re- 
sults for the combined quantities p^ 2 v and p'^v. Usually, the 
pure velocity scaling is considered without density weighting. 
However, for highly supersonic turbulence it is interesting to in- 
vestigate the scaling of combinations of density and velocity. 
Note that CVIs (§ @| and PCA (§ [5]) also analyse convolutions 
of density and velocity statistics. Figure [T| (top panel) shows 
a repetition of Figure [15] (scaling of v) together with the scal- 
ing of p'^ 2 v (middle panel) and p'^v (bottom panel) for direct 
comparison. Since Fourier spectra and A-variance analyses al- 
ways represent the mean squares of these quantities, p'^ 2 v corre- 
sp onds to the sca l ing of the kinetic energy density p v 2 . As show n 
bv lKritsuk et all (120071) fsee also iHenriksenll 1991b iFleckll 19961) . 
p 1/,3 v corresponds to a constant energy flux within the inertial 
range. This idea was first proposed by Lighthiljj (119551) . Using 
the eddy turnover time t( as the typical evolution timescale of a 
turbulent fluctuation on scale I, the constancy of energy flux in 
the inertial range is defined as 

2 2 3 

pv pv pv 

oc oc oc const , (A.l) 

k i/v i 

which leads to the original Kolmo gorovl (I 1 94 ll) scaling (but now 
including density variations), 

p l/ \cc{ 1 ' 3 (A.2) 

for the quantity P 1|f3 v. Us ing the extend ed self- similar i ty hy- 
pothesis (iBenzi et alii 19931) . we showed in lSchmidt et all ([2008) 
that the relative scaling exponents of p'^v provide a more uni- 
versal scaling compared to the velocity scaling without density 
weighting. Figure[T|(bottom panel) shows that the absolute scal- 
ing infOTedJrom theFourier spectra of p^ 3 v is indeed close to 
the [Kolmogorqyl dl94ll) scaling (scaling proportional to k~ 5 l 3 ) 
for solenoid al forcing, which is in agreement with the results 
obtained in Kritsu k et alj d2007t) . However, compressive forc- 
ing yiel ds significantly steeper scaling (also for p ^ 2 v), which is 
close to iBurgersl (I 1 9481) turbulence (scaling proportional to k~ 2 ). 
The corresponding results inferred from the A-variance analyses 
are compatible with the Fourier spectra to within the uncertain- 
ties. Both quantities p^ 2 v and p^ 3 v show breaks in the scaling 
close to the sonic wavenumber k s for compressive forcing. 

Appendix B: Convergence test for the structure 
functions of centroid velocity increments 

For an accurate and reliable determination of the structure func- 
tion scaling, it must be verified that the number of data pairs used 
for sampling the structure functions was high enough to yield 
converged results. There is no general rule to determine a priori 
the number of data pairs necessary, because the required number 
of data pairs depends on the underlying statistics of the mea- 
sured variable itself and on the desired structure function order. 
However, convergence can be tested by increasing the number of 
data pairs used for computing the structure functions. Showing 
that the computed structure functions do not change significantly 
by further increasing the number of data pairs demonstrates con- 
vergence. Furthermore, if convergence is verified for the high- 
est order under consideration, then the structure functions of 
lower order are also converged. This is because the higher-order 
structure functions of a variable q reflect the statistics of higher 
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Fig. .1. Top panels: Same as Figure[l5] Middle panels: Same as top panels, but instead of the Fourier spectra and A-variances of v, 
the Fourier spectra and A-variances of the density-weighted velocity p^ 2 v are shown. The quantity p'' 2 v has physical reference to 
kinetic energy. Bottom panels: Same as middle panels, but the Fourier spectra and A-variances of the density-weighted velocity p^ 3 v 
are shown. The quantity p 1 ^v has physical reference to a constant kinetic energy dissipation within the inertial range (Kritsuk et al. 
2007t ISchmidt et al]l2008h . 



powers of q than the lower order structure functions. This is re- 
flected in the definition of the pth order structure function in 
equation (l26l i. 



Figure iBTI demonstrates convergence for the structure func- 
tions of CVIs with orders p < 6 discussed in § 14.21 We only 
show the compressive forcing case for clarity, but we also ver- 
ified convergence for the solenoidal forcing case with the same 
method. Figure lBTI shows that sampling each structure function 
with roughly 1 .7 x 10 10 data pairs is sufficient to yield converged 
results. The total number of data pairs used to construct the CVI 
structure functions shown in FiguresQT]and[T2]was tnus roughly 
81 x 3 x 1.7 x 10 10 ~ 4.1 x 10 12 from averaging over 81 reali- 
sations of the turbulence and three projections along the x, y and 
z-axes for each of these realisations. 



Appendix C: Resolution study of the Fourier 
spectra and their dependence on the numerical 
scheme 

The resolution and type of numerical method adopted to model 
supersonic turbulence are expected to critically affect the scal- 
ing of Fourier spectrum functions in the inertial rang e (e.g., 
iKlein et all2007l:lKritsuk et alJ2007tlPadoan et alJ2007h . In this 
section, we investigate the dependence of our Fourier spectra on 
the numerical resolution and on the numerical scheme used in 
the present study. 

C.7. Resolution study 

Figure ICTI shows velocity Fourier spectra E(k) defined in equa- 
tion d32l l for numerical resolutions of 256 3 , 512 3 and 1024 3 
grid points. The inertial range scaling is indeed affected by the 
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Fig. B.l. The 1st (p = 1) and 6th (p = 6) order structure func- 
tions of the centroid velocity increments sampled with differ- 
ent numbers of data pairs is shown for a single snapshot at time 
t — 2 T in z-projection for the case of compressive forcing. The 
number of data pairs used for sampling is given in brackets. The 
structure functions of centroid velocity increments are statisti- 
cally converged for p < 6 for sample sizes of at least 1.7 x 10 10 
data pairs per turbulent realisation and per projection as used 
throughout this study. 



one additional run at 512 3 , we switched off the PPM steepening 
algorithm to check its influence on our results. 

Figure IC.2I shows that the velocity spectra E(k) decrease 
faster with increasing diffusion parameter K for wavenumbers 
k > 40. It is expected that the scheme dissipates more ki- 
netic energy on small scales with increasing K, because the 
PPM diffusion algorithm is designed to act on shocks only 
dColella & Woodward! [1984T, eq. 4.5). In contrast, Figure IC21 
demonstrates that the Fourier spectra at wavenumbers k < 40 
are hardly affected by the PPM diffusion algori thm for both 
soleno idal and compressive forcings. Note that Kritsu k et alj 
(2007) reported that their results for the inertial range scaling are 
highly sensitive to the choice of PPM diffusion parameter in the 
ENZO code. However, our results demonstrate that the choice of 
PPM diffusion parameter only affects the inertial range scaling 
within less than 1 %, which is clearly less than the influence of 
the numerical resolution and less than the typical snapshot-to- 
snapshot variations. Figure IC.2I furthermore demonstrates that 
the PPM contact discontinuity steepening has negligible effects 
for simulations of supersonic turbulence. 

The results obtained here support our conclusion in § [6] that 
the Fourier spectra at resolutions of 1024 3 grid cells are robust 
for wavenumbers k < 40. 



numerical resolution. For solenoidal forcing, the inertial range 
scaling exponent j3 at resolution of 256 3 grid cells is roughly 
13% higher than the scaling exponent at a resolution of 1024 3 . 
However, the difference between the inertial range scaling at 
512 3 and 1024 3 is less than 3% for solenoidal forcing. For com- 
pressive forcing, the difference between the inertial range scal- 
ing exponents at resolutions of 512 3 and 1024 3 grid cells is less 
than 1%. This result indicates that the systematic dependence of 
the inertial range scaling on the numerical resolution is less than 
3% for both solenoidal and compressive forcings. It should be 
emphasised that variance effects introduced by different realisa- 
tions of the turbulence are typically on the order of 5-10% (see 
error bars in Figure [T5l). which is higher than the systematic er- 
rors introduced by resolution effects, as long as the numerical 
resolution is at least 512 3 grid cells. 



C.2. Dependence on parameters of the piecewise parabolic 
method 

We used the piecewis e parabolic method (PPM) 
dColella & Woodward! [l984) to integrate the equations of 
hydrodynamics (eqs. [2] and [3]). PPM impro ves on th e finite- 
volume scheme originally developed by Godunov (1959) 
by representing the flow variables with piecewise parabolic 
functions, which makes the PPM second-order accurate in 
smooth flows. However, PPM is also particularly suitable 
for the accurate modelling of turbulent flows involving sharp 
discontinuities, such as shocks and contact discontinuities. For 
that purpose, PPM uses a lower artificial viscosity controlled 
by the PPM diffusion parameter K. In three simulations with 
resolutions of 512 3 grid cells, we varied the PPM diffusion 
parameter K between 0. 0, 0.1 and 0.2. Note that K = 0.1 is the 
value recommended bv lColella & Woodward! (1 19841) . which was 
used for all production runs throughout this study. The PPM 
algorithm furthermore includes a steepening mechanism to keep 
contact discontinuities from spreading over too many cells. In 
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Fig. C.l. Time-averaged velocity Fourier spectra E(k) defined in equation ( l32l for numerical resolutions of 256 3 , 512 3 and 1024 3 
grid points obtained with solenoidal forcing (left) and compressive forcing (right). The inferred inertial range scaling is converged 
to within less than 3% at the typical resolution of 1024 3 grid points used throughout this study for both types of forcing. 




Fig. C.2. Dependence of the time-averaged velocity Fourier spectra E(k) on parameters of the piecewise parabolic method (PPM) 
(IColeila & Woo dward 1984) at fixed resolution of 512 3 grid cells. Varying the PPM diffusion parameter K between 0.0, 0.1 and 0.2 
affects the dissipation range at wavenumbers k > 40. However, the effect of varying the PPM diffusion parameter is negligible for 
k < 40. Switching off the PPM steepening algorithm for contact discontinuities has also virtually no effect on the Fourier spectra at 
k < 40. 



